Here we provide the supplementary material for our study on sexual selection and extinction risk in passerine birds.

Compiling the data

For this project we aimed to compile data on a large set of species in the Passerine bird order that would allow us to account for many evolutionary and environmental predictors of extinction. The following details the source and use of the data.

Speciation Rate

To obtain speciation rate estimates we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) against a pre-existing extinction rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny created by Jetz et al. (2012). We also estimated extinction rate from BAMM.

Sexual Selection

Our proxy for sexual selection was based on measures of sexual dichromatism. Notably, we use two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, which can be divided to obtain a sexual dichromatism index (SDi). The plumage score was scored using the Handbook of the Birds of the World Del Hoyo et al. (n.d.). By scanning images of males and females in this book across multiple patches Dale et al. (2015) were able to obtain mean plumage scores from RGB values. Notably concern about the accuracy of these measurements were raised, thus this dataset was previously compared with a more thorough reflectance dataset of Australian birds, the two datasets correlated but was expectedly noisy. Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008) using colour discriminability as the parameter.

Alongside sexual dichromatism data there is also partial data for other metrics analysed in Dale et al. (2015) such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. We use a phylogenetic PCA (PC1) that estimates male-bias sexual selection from a species level of sexual size dimorphism, social polygyny and paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (through IUCN) (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we were able to obtain estimates of the following (note that not all the data extracted was subsequently used in the analysis):

  • Range Size
  • Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
  • Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
  • Net primary productivity (NPP) estimates and variability
  • Average and standard deviation of human population density in the species ranges.

The code for extracting the environmental data from the species is provided in another R-markdown document Github link, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and se) for each species and exported that data.

#Obtain means
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T )))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(shps.jetz, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) This process was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Generating biologically relevent predictors of extinction

Here we assess the relationship between sexual selection and extinction risk. However, in doing so we attempt to take into account as many other predictors of extinction as possible, primarily through environmental variables. Basically, our model structure seeks to contain:

Extinction/Diversification ~ Sexual selection
+ Range size
+ Short temporal variability of temperature (mean BIOCLIM4)
+ Spatial variability of temperature (PCA1) [residual.PCA1]
+ Long-term variability of temperature (LIG)
+ NPP

These predictors can be obtained from the compiled datasets seen here:

plumage.scores <- read.csv('data/plumage_scores.csv')
#Generate sexual dichromatism score: 
plumage.scores$SDi <- abs(plumage.scores$Male_plumage_score - plumage.scores$Female_plumage_score)

#Make DF with enviro variables and plumage scores
complete.dataframe <- read.csv('data/complete.dataframe.csv')
complete.dataframe <- left_join(plumage.scores %>% dplyr::select(binomial, Male_plumage_score, Female_plumage_score, SDi), complete.dataframe, by = "binomial")

Spatial Variability PCAs

To obtain estimates of spatial variability we can obtain PCAs of the standard errors of the bioclimatic variables (excluding seasonality in temperature and precipitation). From there we can run GAM models with the PCA components and range size and extract residuals. We do this as we expect variability will increase alongside increased range size, thus taking the residuals allows us to correct for this association.

Table S1: Loadings for the first two PCs, when we conduct a PCA on the variation (standard error, where n ~ 1,000) pf each bioclim variable except the standard error of temperature seasonality (se.bioclim4) and the standard error of precipitation seasonality (se.bioclim15), as these are already measures of variability and based on other bioclim variables. We use the standard error as, while we attempted to sample each species range 1,000 times, some points returned NA values and would thus affect the standard deviation, this inconvenience should not affect standard error.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('se.bioclim1', 'se.bioclim2', 'se.bioclim3', 'se.bioclim5', 'se.bioclim6', 'se.bioclim7', 'se.bioclim8', 'se.bioclim9', 'se.bioclim10', 'se.bioclim11', 'se.bioclim12', 'se.bioclim13', 'se.bioclim14', 'se.bioclim16', 'se.bioclim17', 'se.bioclim18', 'se.bioclim19')], #Select se.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()
  PC1 PC2
se.bioclim1 -0.3328 0.08264
se.bioclim2 -0.2383 0.05458
se.bioclim3 -0.2317 -0.06584
se.bioclim5 -0.312 0.08091
se.bioclim6 -0.3288 0.08275
se.bioclim7 -0.2934 0.09748
se.bioclim8 -0.2916 0.1181
se.bioclim9 -0.3187 0.1125
se.bioclim10 -0.3111 0.06704
se.bioclim11 -0.3296 0.08775
se.bioclim12 -0.1386 -0.4246
se.bioclim13 -0.1596 -0.3083
se.bioclim14 -0.01133 -0.3857
se.bioclim16 -0.1708 -0.3254
se.bioclim17 -0.01923 -0.3932
se.bioclim18 -0.1433 -0.3066
se.bioclim19 -0.01797 -0.3815
#run GAM
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-NTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. But not as strong for spatial variability in precipitation (PC2). In the analysis we only used the residuals from PC1

Long term climate variability (LIG anomaly)

To gain estimates for change in climate over the past ~130,000 years we can use the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and prescipitation. Data has also been provided for the last-glacial maximum (LGM) but was not used in the analysis.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absoluted temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()
  PC1.LIG PC2.LIG
bio1.LIG.diff -0.2956 -0.01943
bio2.LIG.diff -0.2442 0.0745
bio3.LIG.diff -0.1436 0.2085
bio5.LIG.diff -0.2867 -0.2566
bio6.LIG.diff -0.4094 0.03346
bio7.LIG.diff -0.4012 -0.09706
bio8.LIG.diff 0.07108 -0.2385
bio9.LIG.diff -0.3365 -0.00753
bio10.LIG.diff -0.07014 -0.3766
bio11.LIG.diff -0.4201 -0.004
bio12.LIG.diff -0.05181 0.3337
bio13.LIG.diff -0.2552 0.2591
bio14.LIG.diff 0.05701 0.3298
bio16.LIG.diff -0.1946 0.3164
bio17.LIG.diff 0.03896 0.3603
bio18.LIG.diff 0.06173 0.277
bio19.LIG.diff 0.08743 0.2856
#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA compontents for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We can check the correlations between environmental predictors that we will use in our model. Specifically we want to know whether the following environmental variables are correlated:

  • Range size
  • Short temporal variability of temperature (mean BIOCLIM4)
  • Spatial variability of temperature (PCA1) [residual.PCA1]
  • Long-term variability of temperature (LIG)
  • NPP

To check some of the correlations between the environmental predictors to be used in the PGLS models we can inspect the following correlation plot with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However this correlation is not excessively high (0.53) so we can be satisfied in including these predictors within our model.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an overdispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interperetation.


Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

Diversification rates can be simply estimated from thew ‘tippiness of the phylogeny’. Here a log equal splits (also called DR) value can be achieved that is weighted to more recent diversification. It was used by Jetz et al. (2012) and recently (in comparison with BAMM) in a study ivestigation bird diversity alongside elevation (Quintero and Jetz 2018). We can look at speciation rate using ES-Sim adapted by Harvey Michael et al. (2017) and found here https://github.com/mgharvey/ES-sim.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 100 trees plus a full MCC tree of passerine birds based on 2500 samples of the posterior (Jetz et al. 2012). This was done using the phanghorn (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('bamm_extinction/data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

To gain an estimate of variation between trees we can take the values and obtain a mean, variance and coeeficient of variation (CV). We obtain a CV value to account for unequal variation at higher estimates. The CV is expressed as the ratio of the standard deviation (\(\sigma\)) to the mean (\(\mu\)) (converting to a percentage not needed):

\(C_V = \frac{\sigma}{\mu} \times 100\)

Furthermore, we can account for overdispersion by obtaining a transformed version of CV. In such a case we can calculate CV for a log-normal distribution:

\({\widehat {c_{\rm {v}}}}_{\rm {ln}}={\sqrt {\mathrm {e} ^{s_{\rm {ln}}^{2}}-1}}\)

Given that our ES value is log-transformed, the calculation of \({\widehat {c_{\rm {v}}}}_{\rm {ln}}\) is perhaps more ideal. As for ND we can see from the a qqnorm plot that value of ND is normally distributed, therefore the regular CV is more appropriate for this distribution.

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')
dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = mean.loges, y = mean.nd, 
                          fill = (1/(CV.nd))/sum((1/(CV.nd)))*100, 
                      size = (1/(CV.loges))/sum((1/(CV.loges)))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.loges - 0.674*sqrt(var.loges), xmax = mean.loges + 0.674*sqrt(var.loges)), 
                 size = 0.015)+
  geom_errorbar(aes(ymin = mean.nd - 0.674*sqrt(var.nd), ymax = mean.nd + 0.674*sqrt(var.nd)), 
                size = 0.015)+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using DR]', y='Node Density [ND]', x= 'Log Diversification Rate [DR]', fill = 'Weight (%) [using ND] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1, breaks=c(0,0.015,.03),
                           limits=c(0,0.03))
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)
grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) are correlated in their estimation of species rate. Furthermore, we see similar variation across each species in the estimates across 100 trees as the inverse of the coefficeint of variation (weight) is approximately equivalent across all values of ND and DR (colours and size). Notably, these estimates speciation rate are relatively uncertain across all species, with 50 % CIs in both axes plotted.

Combine the estimates of DR and ND with the restricted dataframe containing all the other variables. We can inspect the trend in SDi in relation to the speciation rates.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures the pattern and spread is similar, with no obvious relationship. Note that terminal branch lenght (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

Set up BAMM parameters

For the use of BAMM we used the following code to generate the parameters across the 100 trees. Each paramaeter value is specified in this code chunk to ensure reproducibility when re-running on any of the trees. These same parameters were also used for the MCC tree (with the seed set at 2500).

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

Run analysis

BAMM can be run through the terminal through the following syntax: bamm -c control_tree_xxxx.txt. To generate these commands we can use a loop function, from which we get:

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

The analysis was run over multiple CPU’s, each generating a respective MCMC and EventData output. Due to the size of the event data file (~ 50 Gb in total) they are not included as supplementary material here. However we have simplified the event data objects into tip rate estimates of the mean and variance across 100 trees + MCC.

Read in the event data and extract tip data

Now we have a series of trees and event data we can read in. Firstly and most importantly let’s check for convergence in the MCC tree through the following:

Table S3: Effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) for the run on the MCC indicate that BAMM converges (ESS > 200).

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()
  N_Shifts logLik
Effective Sample Size 2456 1302

Those values are well above 200, so we are good to use the MCC data in analysis! We can also check the convergence of BAMM across 100 runs of BAMM. The Raw MCMCs are not included in this file or repository but we can read in a dataframe that has extracted the ESS for each of the runs.

Table S4: For the 100 trees that BAMM was run on effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) also indicates that BAMM converges with the minimum for each parameter across the 100 trees being over 200.

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()
N_Shifts logLik
Min. : 723.4 Min. : 278.7
1st Qu.:1396.1 1st Qu.: 511.3
Median :1618.6 Median : 638.6
Mean :1709.5 Mean : 669.9
3rd Qu.:2070.8 3rd Qu.: 804.8
Max. :3425.6 Max. :1244.1

Given that it appears BAMM converges we can we can make a data frame with the mean and variance for extinction and speciation tip rates from the large event data set for the MCC with the following code and then plot the variation we see in tip-rates.

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)
## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence
saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(x)
  var.lambda = var(x)
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(x)
  var.mu = var(x)
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = log(mean.lambda), y = log(mean.mu), 
                          fill = (1/(log(var.lambda)))/sum((1/(log(var.lambda))))*100, 
                      size = (1/(log(var.mu)))/sum((1/(log(var.mu))))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = log(mean.lambda) - 0.674*log(sqrt(var.lambda)), xmax = log(mean.lambda) + 0.674*log(sqrt(var.lambda))), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = log(mean.mu) - 0.674*log(sqrt(var.mu)), ymax = log(mean.mu) + 0.674*log(sqrt(var.mu))), 
                size = 0.0025)+
  xlim(-8, 3)+
  ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight (%) [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1, breaks=c(0,.15,0.3),
                           limits=c(-0,0.3))

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)
grid.newpage()
grid.draw(BAMM.variance)

Figure S5: The tip-rate estimate for BAMM are highly variable within each run of BAMM. Across most species there is high variability in the posterior distribution of tip-rate estimates. Here we show mean values, weights based on the variance and 50 % CIs.

Analysis of BAMM results

Given the high variability in tip-rate estimates from the above plot, below we peform some diagnostics on BAMM to demonstrate that the variability is unlikely an error in sampling or parameters, rather an inherent aspect of BAMM. The following is not inherently necessary to understandfor the conclusions drawn from the paper, however they do raise a set of methodological questions about BAMM that warrant further investigation.

Credible number of shifts

To plot the credible shift set, we need the prior distribution on the number of rate shifts (this is generated internally by BAMMtools). We can then estimate the credible set of rate shifts using the BAMMtools function credibleShiftSet:

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css$number.distinct
## [1] 950

This value indicates that from a posterior sample of 1000 we have 950 unique rate shifts, which is very high. Arguably, across such a large phylogenetic tree, even if the number of shifts reaches convergence the location of those > 50 shifts will be quite variable, with the variation likely impacting tip-rate estimates. Even if we increase the sampling of the posterior to 2,500 instead of 1,000 we still see that the CSS is high (close to 2500). We can have a look what our top few shifts are:

summary(css)
## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

Notably, each of the top shifts has low probability, indicating that out of the entire posterior sample we cannot differentiate which shift configuration is more likely. Based on the follwing eastimates of the BayesFactor, we are convinced that the number of shifts is non-zero, however it is still quite a wide distribution, sourced from potentially weakly non-informative priors. But, this is at odds with the ESS and the following plot:

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)
##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56
plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

Figure S6: The apparant convergence of the number of shifts in the posterior is at odds with the variability seen in the CSS and although there is greater certainty in the number of shifts being within the range above, the position of those shifts remain highly variable.

We can compare our run of BAMM against Harvey et al. (2017) who used BAMM on a genetic-only MCC tree with different parameters

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)
css3$number.distinct
## [1] 1188

Based on Their posterior sample of 1250, 1188 are unique. Like our run on the MCC, we find high uncertainity in the position of shift-configurations.

Let’s run a model based on genetic data from the Harvey MCC

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = log(mean.lambda), y = log(mean.mu), 
                          fill = (1/(log(var.lambda)))/sum((1/(log(var.lambda))))*100, 
                      size = (1/(log(var.mu)))/sum((1/(log(var.mu))))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = log(mean.lambda) - 0.674*log(sqrt(var.lambda)), xmax = log(mean.lambda) + 0.674*log(sqrt(var.lambda))), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = log(mean.mu) - 0.674*log(sqrt(var.mu)), ymax = log(mean.mu) + 0.674*log(sqrt(var.mu))), 
                size = 0.0025)+
  xlim(-8, 3)+
  ylim(-12.5,1)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight (%) [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1, breaks=c(0,.15,0.3),
                           limits=c(-0,0.3))

# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

Figure S7: Based on the event data from the BAMM run by Harvey et al. (2017) we find that in both our case and theirs the tip-rate estimates for many species are extrememly variable. here we present mean estimates with 50 % CIs.

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive. Based on preliminary analysis even slight differences in lambda lead to large reductions in AIC.

2: Create a global model of six predictor variables plus the fiver interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

PGLS Models on DR and ND

Using corPagel we can estimate the phylogenetic signal for a model with all predictors (interactions do not appear to affect the estimate of \(\lambda\)):

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

Inspect the \(\lambda\) value, which can then be fixed for successive models.

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.985
MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()
  • ND lambda: 0.9996

The lambda is high and not too dissimilar than if we assume Brownian Motion \(\lambda = 1\). However, given the large sample size this difference may have an effect on the results so, to play safe we included it as a fixed value for \(\lambda\) in all successive models.

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

Table S5: The dredged models both show the top model is one with no interactions, with \(\delta AICc > 4\) in both cases. We can be resonably confident that interactions are unlikely to affect the pattern of speciation we see in passerine birds.

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.879811 2.14e-05 -0.0065780 3e-07 0.0015864 -0.0024295 -0.0012787 NA NA NA NA NA 8 -3409.547 6835.118 0.00000 9.988652e-01
8 -2.879495 2.14e-05 -0.0065788 3e-07 -0.0001028 -0.0023765 -0.0012942 NA NA NA 0.0003487 NA 9 -3416.089 6850.208 15.09006 5.281325e-04
16 -2.879345 2.15e-05 -0.0066034 3e-07 0.0015509 -0.0015982 -0.0013082 NA NA NA NA -0.0001638 9 -3416.488 6851.007 15.88895 3.542138e-04
2 -2.884345 2.14e-05 -0.0064052 3e-07 0.0015991 -0.0024487 -0.0003035 NA -0.0000369 NA NA NA 9 -3416.842 6851.714 16.59633 2.486904e-04
1 -2.878432 1.57e-05 -0.0065688 3e-07 0.0015872 -0.0024180 -0.0015875 1.2e-06 NA NA NA NA 9 -3421.210 6860.452 25.33361 3.150533e-06
24 -2.878703 2.16e-05 -0.0066190 3e-07 -0.0004658 -0.0010549 -0.0013437 NA NA NA 0.0004121 -0.0002586 10 -3422.848 6865.734 30.61631 2.245227e-07
4 -2.884812 2.14e-05 -0.0065580 7e-07 0.0015560 -0.0024009 -0.0003587 NA NA -1e-07 NA NA 9 -3424.307 6866.644 31.52618 1.424567e-07
10 -2.890455 2.13e-05 -0.0061602 3e-07 -0.0001797 -0.0024196 0.0010665 NA -0.0000894 NA 0.0003710 NA 10 -3423.322 6866.681 31.56307 1.398533e-07
18 -2.885405 2.15e-05 -0.0063727 3e-07 0.0015667 -0.0015933 -0.0000022 NA -0.0000495 NA NA -0.0001699 10 -3423.771 6867.580 32.46202 8.922136e-08
9 -2.879477 2.13e-05 -0.0065786 3e-07 -0.0000996 -0.0023764 -0.0012983 0.0e+00 NA NA 0.0003480 NA 10 -3427.745 6875.529 40.41070 1.676618e-09
17 -2.876970 1.25e-05 -0.0065992 3e-07 0.0015377 -0.0012411 -0.0018098 1.9e-06 NA NA NA -0.0002306 10 -3428.008 6876.054 40.93602 1.289324e-09
3 -2.889434 1.34e-05 -0.0061257 3e-07 0.0016199 -0.0024624 0.0007743 1.7e-06 -0.0000939 NA NA NA 10 -3428.372 6876.782 41.66339 8.962234e-10
12 -2.884227 2.14e-05 -0.0065599 7e-07 -0.0000916 -0.0023507 -0.0004248 NA NA -1e-07 0.0003405 NA 10 -3430.871 6881.780 46.66167 7.362996e-11
26 -2.893478 2.15e-05 -0.0060556 3e-07 -0.0006022 -0.0009958 0.0018498 NA -0.0001211 NA 0.0004478 -0.0002816 11 -3430.028 6882.101 46.98312 6.269768e-11
20 -2.884117 2.15e-05 -0.0065817 7e-07 0.0015263 -0.0016648 -0.0004369 NA NA -1e-07 NA -0.0001454 10 -3431.269 6882.576 47.45753 4.945782e-11
6 -2.895008 2.14e-05 -0.0061863 7e-07 0.0015803 -0.0024391 0.0018191 NA -0.0000790 -1e-07 NA NA 10 -3431.545 6883.128 48.01036 3.751366e-11
25 -2.877827 1.81e-05 -0.0066164 3e-07 -0.0003444 -0.0009513 -0.0015348 7.0e-07 NA NA 0.0003863 -0.0002784 11 -3434.448 6890.942 55.82373 7.542887e-13
11 -2.891727 1.89e-05 -0.0060861 3e-07 -0.0000889 -0.0024252 0.0013324 5.0e-07 -0.0001044 NA 0.0003533 NA 11 -3434.894 6891.834 56.71575 4.828789e-13
5 -2.883424 1.60e-05 -0.0065495 7e-07 0.0015572 -0.0023903 -0.0006670 1.1e-06 NA -1e-07 NA NA 10 -3435.976 6891.989 56.87126 4.467550e-13
19 -2.894722 8.10e-06 -0.0058780 3e-07 0.0015805 -0.0010549 0.0020214 2.8e-06 -0.0001542 NA NA -0.0002814 11 -3435.056 6892.158 57.03960 4.106919e-13
28 -2.882980 2.16e-05 -0.0065992 6e-07 -0.0004288 -0.0011303 -0.0005650 NA NA -1e-07 0.0004000 -0.0002393 11 -3437.666 6897.377 62.25931 3.020435e-14
14 -2.901190 2.13e-05 -0.0059396 7e-07 -0.0002030 -0.0024099 0.0032037 NA -0.0001318 -1e-07 0.0003719 NA 11 -3438.023 6898.091 62.97260 2.114363e-14
22 -2.895506 2.14e-05 -0.0061663 7e-07 0.0015517 -0.0016637 0.0019999 NA -0.0000885 -1e-07 NA -0.0001541 11 -3438.494 6899.034 63.91561 1.319496e-14
27 -2.898432 1.32e-05 -0.0057816 3e-07 -0.0003798 -0.0007224 0.0029159 1.7e-06 -0.0001787 NA 0.0004037 -0.0003394 12 -3441.451 6906.955 71.83674 2.513946e-16
13 -2.884235 2.15e-05 -0.0065599 7e-07 -0.0000930 -0.0023508 -0.0004230 0.0e+00 NA -1e-07 0.0003408 NA 11 -3442.528 6907.101 71.98280 2.336900e-16
21 -2.881649 1.31e-05 -0.0065788 6e-07 0.0015150 -0.0013253 -0.0009513 1.8e-06 NA -1e-07 NA -0.0002092 11 -3442.806 6907.657 72.53841 1.770074e-16
7 -2.901278 1.25e-05 -0.0058625 8e-07 0.0016024 -0.0024539 0.0031392 1.8e-06 -0.0001449 -1e-07 NA NA 11 -3443.046 6908.138 73.01974 1.391463e-16
30 -2.903249 2.15e-05 -0.0058572 7e-07 -0.0005994 -0.0010691 0.0037845 NA -0.0001586 -1e-07 0.0004443 -0.0002653 12 -3444.762 6913.578 78.46031 9.163629e-18
29 -2.882138 1.85e-05 -0.0065971 6e-07 -0.0003202 -0.0010362 -0.0007469 7.0e-07 NA -1e-07 0.0003770 -0.0002574 12 -3449.269 6922.592 87.47390 1.011093e-19
15 -2.903204 1.79e-05 -0.0058308 7e-07 -0.0000759 -0.0024175 0.0036222 7.0e-07 -0.0001538 -1e-07 0.0003472 NA 12 -3449.583 6923.220 88.10200 7.385839e-20
23 -2.905869 7.40e-06 -0.0056354 7e-07 0.0015654 -0.0011030 0.0042394 2.9e-06 -0.0002007 -1e-07 NA -0.0002702 12 -3449.752 6923.557 88.43915 6.240085e-20
31 -2.909040 1.25e-05 -0.0055511 7e-07 -0.0003574 -0.0007748 0.0050239 1.9e-06 -0.0002227 -1e-07 0.0003961 -0.0003275 13 -3456.166 6938.395 103.27731 3.742176e-23
kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 0.1404205 5e-07 -0.0001462 0 -0.0000102 -6.47e-05 -0.0000575 NA NA NA NA NA 8 14329.64 -28643.26 0.00000 9.999550e-01
8 0.1403564 4e-07 -0.0001422 0 -0.0000850 -5.71e-05 -0.0000618 NA NA NA 1.48e-05 NA 9 14320.05 -28622.06 21.19677 2.495513e-05
16 0.1404275 5e-07 -0.0001466 0 -0.0000108 -5.67e-05 -0.0000581 NA NA NA NA -1.6e-06 9 14319.25 -28620.46 22.79683 1.121276e-05
2 0.1400429 5e-07 -0.0001322 0 -0.0000107 -6.61e-05 0.0000205 NA -2.9e-06 NA NA NA 9 14318.99 -28619.95 23.30593 8.692847e-06
1 0.1403637 6e-07 -0.0001457 0 -0.0000101 -6.47e-05 -0.0000496 0e+00 NA NA NA NA 9 14314.58 -28611.13 32.12698 1.056075e-07
4 0.1406803 5e-07 -0.0001499 0 -0.0000090 -6.52e-05 -0.0000932 NA NA 0 NA NA 9 14311.61 -28605.19 38.06105 5.434109e-09
24 0.1403735 4e-07 -0.0001429 0 -0.0000910 -3.30e-05 -0.0000638 NA NA NA 1.57e-05 -4.8e-06 10 14309.74 -28599.43 43.82308 3.047333e-10
10 0.1398108 4e-07 -0.0001218 0 -0.0000886 -5.88e-05 0.0000503 NA -4.2e-06 NA 1.54e-05 NA 10 14309.46 -28598.88 44.37966 2.307068e-10
18 0.1400302 5e-07 -0.0001318 0 -0.0000116 -5.60e-05 0.0000242 NA -3.1e-06 NA NA -2.1e-06 10 14308.61 -28597.18 46.08010 9.858583e-11
9 0.1401845 8e-07 -0.0001400 0 -0.0000977 -5.59e-05 -0.0000401 -1e-07 NA NA 1.74e-05 NA 10 14305.25 -28590.46 52.79550 3.432282e-12
17 0.1403708 6e-07 -0.0001459 0 -0.0000104 -6.13e-05 -0.0000504 0e+00 NA NA NA -7.0e-07 10 14304.23 -28588.42 54.83371 1.238774e-12
3 0.1400821 5e-07 -0.0001346 0 -0.0000106 -6.59e-05 0.0000102 0e+00 -2.4e-06 NA NA NA 10 14304.00 -28587.96 55.29503 9.835981e-13
12 0.1406450 4e-07 -0.0001461 0 -0.0000867 -5.74e-05 -0.0001020 NA NA 0 1.54e-05 NA 10 14302.09 -28584.14 59.11661 1.455366e-13
20 0.1406914 5e-07 -0.0001503 0 -0.0000097 -5.54e-05 -0.0000943 NA NA 0 NA -2.0e-06 10 14301.22 -28582.41 60.84729 6.125736e-14
6 0.1404437 5e-07 -0.0001414 0 -0.0000094 -6.61e-05 -0.0000452 NA -1.7e-06 0 NA NA 10 14300.95 -28581.85 61.40235 4.641190e-14
26 0.1397605 4e-07 -0.0001201 0 -0.0000961 -3.08e-05 0.0000624 NA -4.8e-06 NA 1.64e-05 -5.7e-06 11 14299.18 -28576.32 66.93341 2.921291e-15
5 0.1406261 6e-07 -0.0001493 0 -0.0000090 -6.53e-05 -0.0000857 0e+00 NA 0 NA NA 10 14296.55 -28573.06 70.19982 5.705361e-16
25 0.1402082 8e-07 -0.0001405 0 -0.0000999 -4.30e-05 -0.0000430 -1e-07 NA NA 1.76e-05 -2.6e-06 11 14294.92 -28567.80 75.45446 4.123369e-17
11 0.1399495 7e-07 -0.0001307 0 -0.0000978 -5.69e-05 0.0000099 -1e-07 -2.0e-06 NA 1.73e-05 NA 11 14294.66 -28567.27 75.98128 3.168506e-17
19 0.1400559 5e-07 -0.0001333 0 -0.0000113 -5.77e-05 0.0000174 0e+00 -2.7e-06 NA NA -1.7e-06 11 14293.70 -28565.35 77.90868 1.208721e-17
28 0.1406725 4e-07 -0.0001470 0 -0.0000935 -3.05e-05 -0.0001055 NA NA 0 1.64e-05 -5.4e-06 11 14291.80 -28561.55 81.70210 1.813827e-18
14 0.1402430 4e-07 -0.0001316 0 -0.0000890 -5.86e-05 -0.0000209 NA -3.0e-06 0 1.58e-05 NA 11 14291.46 -28560.87 82.38244 1.290804e-18
22 0.1404322 5e-07 -0.0001410 0 -0.0000103 -5.50e-05 -0.0000416 NA -1.9e-06 0 NA -2.3e-06 11 14290.57 -28559.08 84.17121 5.277578e-19
13 0.1404724 8e-07 -0.0001438 0 -0.0000991 -5.62e-05 -0.0000802 -1e-07 NA 0 1.79e-05 NA 11 14287.28 -28552.52 90.73241 1.984664e-20
21 0.1406404 6e-07 -0.0001497 0 -0.0000094 -5.93e-05 -0.0000874 0e+00 NA 0 NA -1.2e-06 11 14286.20 -28550.36 92.89788 6.721422e-21
7 0.1404979 6e-07 -0.0001445 0 -0.0000092 -6.58e-05 -0.0000590 0e+00 -1.0e-06 0 NA NA 11 14285.96 -28549.88 93.37861 5.285327e-21
27 0.1398914 7e-07 -0.0001279 0 -0.0001008 -3.94e-05 0.0000252 -1e-07 -2.8e-06 NA 1.76e-05 -3.6e-06 12 14284.39 -28544.73 98.52951 4.023185e-22
30 0.1401987 4e-07 -0.0001299 0 -0.0000970 -2.90e-05 -0.0000096 NA -3.5e-06 0 1.69e-05 -6.0e-06 12 14281.20 -28538.34 104.91315 1.653402e-23
29 0.1405076 8e-07 -0.0001446 0 -0.0001019 -4.01e-05 -0.0000847 -1e-07 NA 0 1.82e-05 -3.3e-06 12 14276.97 -28529.89 113.36490 2.416050e-25
15 0.1404182 8e-07 -0.0001418 0 -0.0000991 -5.65e-05 -0.0000689 -1e-07 -4.0e-07 0 1.79e-05 NA 12 14276.70 -28529.34 113.91919 1.831227e-25
23 0.1404716 5e-07 -0.0001432 0 -0.0000099 -5.75e-05 -0.0000518 0e+00 -1.4e-06 0 NA -1.7e-06 12 14275.66 -28527.26 115.99255 6.494083e-26
31 0.1403601 7e-07 -0.0001390 0 -0.0001023 -3.86e-05 -0.0000535 -1e-07 -1.2e-06 0 1.82e-05 -3.7e-06 13 14266.43 -28506.79 136.46456 2.328506e-30

From the dredged model list we can see that the top model is one with no interactions. Run this model on the MCC and 100 trees. Noting that each model uses a unique set of values for DR/ND and a unique tree in the correlation structure.

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
DR.model.data <- lapply(DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
DR.pgls.models <- mcmapply(function(x,y) {
  gls(ES ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.985, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = DR.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")

#Now for Node Density:
ND.model.data <- lapply(nd.list, function(x) {
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
ND.model.data <- lapply(ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.9996, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = ND.model.data, y = pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

From these models we can get a distribution of the estimates for the models alongside the 95 % CIs of the MCC model, this will show the variation beween trees, and the variation associated with the top MCC model.

DR.pgls.models <- readRDS("data/DR.pgls.models.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')

DR.pgls.summary <- lapply(DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.DR.summary <- data.frame(MCC.DR.top$coefficients, confint(MCC.DR.top)) %>% tibble::rownames_to_column()

DR.pgls.summary <- bind_rows(DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")




#Now plot this for Node Density (ND):  

ND.pgls.models <- readRDS("data/ND.pgls.models.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')

ND.pgls.summary <- lapply(ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.ND.summary <- data.frame(MCC.ND.top$coefficients, confint(MCC.ND.top)) %>% tibble::rownames_to_column()

ND.pgls.summary <- bind_rows(ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we can now undertake model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, meaning more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagal’s lambda was = 1 and running corPagel lead to problemd of convergence. Therefore the following models are run assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(log(mean.lambda) ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(log(mean.mu) ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Speciation Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.648503 -2.0e-06 -1.37e-05 0 -1.12e-05 3.50e-05 0.0000584 NA NA NA NA NA 8 8353.254 -16690.48 0.00000 9.999123e-01
8 -2.648548 -2.0e-06 -1.18e-05 0 -6.30e-05 4.19e-05 0.0000554 NA NA NA 1.13e-05 NA 9 8343.888 -16669.74 20.73974 3.136067e-05
16 -2.648489 -2.0e-06 -1.40e-05 0 -1.15e-05 4.98e-05 0.0000573 NA NA NA NA -3.0e-06 9 8343.836 -16669.64 20.84310 2.978112e-05
2 -2.647505 -2.0e-06 -5.18e-05 0 -8.30e-06 3.69e-05 -0.0001854 NA 9.20e-06 NA NA NA 9 8343.708 -16669.39 21.09789 2.621880e-05
1 -2.648473 -2.0e-06 -1.42e-05 0 -1.18e-05 3.51e-05 0.0000548 0e+00 NA NA NA NA 9 8339.205 -16660.38 30.10536 2.901787e-07
4 -2.648288 -2.0e-06 -1.68e-05 0 -1.11e-05 3.85e-05 0.0000279 NA NA 0 NA NA 9 8336.072 -16654.11 36.36989 1.265731e-08
24 -2.648526 -2.0e-06 -1.20e-05 0 -7.13e-05 7.46e-05 0.0000525 NA NA NA 1.30e-05 -6.5e-06 10 8334.520 -16649.00 41.48242 9.821319e-10
10 -2.647644 -2.1e-06 -4.63e-05 0 -5.57e-05 4.30e-05 -0.0001645 NA 8.30e-06 NA 1.03e-05 NA 10 8334.336 -16648.63 41.85004 8.172261e-10
18 -2.647515 -2.0e-06 -5.12e-05 0 -8.50e-06 4.60e-05 -0.0001817 NA 9.00e-06 NA NA -1.9e-06 10 8334.293 -16648.55 41.93551 7.830375e-10
9 -2.648609 -1.9e-06 -1.06e-05 0 -6.59e-05 4.24e-05 0.0000621 0e+00 NA NA 1.22e-05 NA 10 8329.896 -16639.75 50.72978 9.641194e-12
17 -2.648432 -2.1e-06 -1.50e-05 0 -1.25e-05 5.34e-05 0.0000505 0e+00 NA NA NA -3.8e-06 10 8329.834 -16639.63 50.85452 9.058213e-12
3 -2.647431 -1.8e-06 -5.65e-05 0 -6.20e-06 3.71e-05 -0.0002150 0e+00 1.07e-05 NA NA NA 10 8329.770 -16639.50 50.98172 8.500083e-12
12 -2.648321 -2.0e-06 -1.50e-05 0 -6.41e-05 4.57e-05 0.0000228 NA NA 0 1.16e-05 NA 10 8326.709 -16633.38 57.10299 3.982947e-13
20 -2.648279 -2.0e-06 -1.70e-05 0 -1.13e-05 5.17e-05 0.0000274 NA NA 0 NA -2.7e-06 10 8326.654 -16633.27 57.21322 3.769381e-13
6 -2.647016 -2.0e-06 -6.33e-05 0 -7.50e-06 4.20e-05 -0.0002741 NA 1.09e-05 0 NA NA 10 8326.570 -16633.10 57.38149 3.465218e-13
26 -2.647690 -2.1e-06 -4.41e-05 0 -6.28e-05 6.88e-05 -0.0001512 NA 7.70e-06 NA 1.17e-05 -5.2e-06 11 8324.971 -16627.90 62.58783 2.565593e-14
5 -2.648251 -2.0e-06 -1.75e-05 0 -1.17e-05 3.85e-05 0.0000233 0e+00 NA 0 NA NA 10 8322.024 -16624.01 66.47377 3.675915e-15
25 -2.648553 -2.0e-06 -1.15e-05 0 -7.23e-05 7.34e-05 0.0000556 0e+00 NA NA 1.33e-05 -6.2e-06 11 8320.553 -16619.06 71.42387 3.093612e-16
11 -2.647523 -1.7e-06 -5.46e-05 0 -6.27e-05 4.48e-05 -0.0002202 -1e-07 1.12e-05 NA 1.28e-05 NA 11 8320.468 -16618.89 71.59252 2.843436e-16
19 -2.647437 -1.8e-06 -5.61e-05 0 -6.40e-06 3.94e-05 -0.0002125 0e+00 1.06e-05 NA NA -5.0e-07 11 8320.432 -16618.82 71.66592 2.740978e-16
28 -2.648306 -2.0e-06 -1.52e-05 0 -7.20e-05 7.69e-05 0.0000211 NA NA 0 1.32e-05 -6.2e-06 11 8317.341 -16612.64 77.84794 1.245936e-17
14 -2.647149 -2.1e-06 -5.80e-05 0 -5.56e-05 4.83e-05 -0.0002543 NA 1.00e-05 0 1.05e-05 NA 11 8317.199 -16612.35 78.13101 1.081507e-17
22 -2.647027 -2.0e-06 -6.29e-05 0 -7.60e-06 4.80e-05 -0.0002708 NA 1.08e-05 0 NA -1.2e-06 11 8317.156 -16612.27 78.21677 1.036109e-17
13 -2.648379 -1.9e-06 -1.39e-05 0 -6.67e-05 4.61e-05 0.0000293 0e+00 NA 0 1.24e-05 NA 11 8312.718 -16603.39 87.09376 1.224030e-19
21 -2.648217 -2.1e-06 -1.82e-05 0 -1.23e-05 5.56e-05 0.0000199 0e+00 NA 0 NA -3.5e-06 11 8312.653 -16603.26 87.22383 1.146959e-19
7 -2.646914 -1.8e-06 -6.91e-05 0 -5.00e-06 4.24e-05 -0.0003112 0e+00 1.27e-05 0 NA NA 11 8312.636 -16603.23 87.25674 1.128240e-19
27 -2.647565 -1.8e-06 -5.23e-05 0 -6.59e-05 5.94e-05 -0.0002056 -1e-07 1.05e-05 NA 1.33e-05 -3.0e-06 12 8311.150 -16598.25 92.23739 9.351196e-21
30 -2.647204 -2.1e-06 -5.56e-05 0 -6.18e-05 7.07e-05 -0.0002400 NA 9.50e-06 0 1.17e-05 -4.5e-06 12 8307.834 -16591.61 98.87041 3.392537e-22
29 -2.648331 -2.0e-06 -1.47e-05 0 -7.28e-05 7.60e-05 0.0000239 0e+00 NA 0 1.35e-05 -6.0e-06 12 8303.374 -16582.69 107.78986 3.923668e-24
15 -2.646985 -1.7e-06 -6.78e-05 0 -6.32e-05 5.06e-05 -0.0003213 -1e-07 1.34e-05 0 1.32e-05 NA 12 8303.340 -16582.63 107.85783 3.792547e-24
23 -2.646904 -1.8e-06 -6.96e-05 0 -4.80e-06 3.96e-05 -0.0003147 -1e-07 1.29e-05 0 NA 6.0e-07 12 8303.303 -16582.55 107.93102 3.656279e-24
31 -2.647020 -1.7e-06 -6.61e-05 0 -6.53e-05 6.00e-05 -0.0003103 -1e-07 1.29e-05 0 1.35e-05 -1.9e-06 13 8294.024 -16561.99 128.49799 1.250196e-28
kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Extinction Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -6.812225 4e-07 -0.0000888 0 -0.0001928 -1.94e-05 0.0000361 NA NA NA NA NA 8 5717.114 -11418.20 0.00000 9.998665e-01
16 -6.812223 4e-07 -0.0000891 0 -0.0001934 2.40e-06 0.0000343 NA NA NA NA -4.7e-06 9 5708.174 -11398.32 19.88684 4.803637e-05
8 -6.812229 4e-07 -0.0000891 0 -0.0001779 -2.15e-05 0.0000371 NA NA NA -3.4e-06 NA 9 5708.136 -11398.24 19.96336 4.623318e-05
2 -6.811955 4e-07 -0.0000991 0 -0.0001921 -1.88e-05 -0.0000324 NA 2.6e-06 NA NA NA 9 5707.959 -11397.89 20.31734 3.873367e-05
1 -6.812299 5e-07 -0.0000876 0 -0.0001918 -1.94e-05 0.0000442 0e+00 NA NA NA NA 9 5703.513 -11388.99 29.20883 4.542818e-07
4 -6.812053 4e-07 -0.0000912 0 -0.0001928 -1.74e-05 0.0000132 NA NA 0 NA NA 9 5700.350 -11382.67 35.53393 1.922416e-08
24 -6.812226 4e-07 -0.0000893 0 -0.0001827 -2.00e-06 0.0000352 NA NA NA -2.4e-06 -4.1e-06 10 5699.228 -11378.42 39.78584 2.293804e-09
18 -6.811993 4e-07 -0.0000979 0 -0.0001928 1.60e-06 -0.0000241 NA 2.2e-06 NA NA -4.4e-06 10 5699.024 -11378.01 40.19459 1.869811e-09
10 -6.811924 4e-07 -0.0001008 0 -0.0001755 -2.11e-05 -0.0000405 NA 2.9e-06 NA -3.8e-06 NA 10 5698.987 -11377.94 40.26855 1.801929e-09
17 -6.812267 5e-07 -0.0000884 0 -0.0001927 -2.00e-07 0.0000394 0e+00 NA NA NA -4.1e-06 10 5694.618 -11369.20 49.00476 2.283987e-11
9 -6.812280 5e-07 -0.0000882 0 -0.0001801 -2.11e-05 0.0000425 0e+00 NA NA -2.7e-06 NA 10 5694.589 -11369.14 49.06300 2.218440e-11
3 -6.811882 6e-07 -0.0001049 0 -0.0001899 -1.84e-05 -0.0000705 -1e-07 4.5e-06 NA NA NA 10 5694.467 -11368.90 49.30697 1.963677e-11
20 -6.812053 4e-07 -0.0000915 0 -0.0001933 4.00e-06 0.0000117 NA NA 0 NA -4.6e-06 10 5691.410 -11362.78 55.42135 9.233120e-13
12 -6.812060 4e-07 -0.0000914 0 -0.0001784 -1.95e-05 0.0000144 NA NA 0 -3.3e-06 NA 10 5691.372 -11362.71 55.49754 8.888011e-13
6 -6.811629 3e-07 -0.0001066 0 -0.0001918 -1.62e-05 -0.0000908 NA 3.7e-06 0 NA NA 10 5691.222 -11362.41 55.79688 7.652497e-13
26 -6.811963 4e-07 -0.0000994 0 -0.0001801 -3.80e-06 -0.0000317 NA 2.5e-06 NA -2.8e-06 -3.6e-06 11 5690.088 -11358.13 60.07356 9.018520e-14
5 -6.812126 5e-07 -0.0000901 0 -0.0001918 -1.75e-05 0.0000211 0e+00 NA 0 NA NA 10 5686.750 -11353.46 64.74224 8.736630e-15
25 -6.812256 4e-07 -0.0000887 0 -0.0001837 -3.20e-06 0.0000385 0e+00 NA NA -2.1e-06 -3.8e-06 11 5685.710 -11349.37 68.82949 1.131900e-15
19 -6.811923 5e-07 -0.0001028 0 -0.0001908 -5.10e-06 -0.0000565 0e+00 3.9e-06 NA NA -2.9e-06 11 5685.611 -11349.18 69.02734 1.025288e-15
11 -6.811873 6e-07 -0.0001051 0 -0.0001791 -2.00e-05 -0.0000698 0e+00 4.5e-06 NA -2.5e-06 NA 11 5685.544 -11349.04 69.16035 9.593193e-16
28 -6.812058 4e-07 -0.0000916 0 -0.0001832 -3.00e-07 0.0000127 NA NA 0 -2.3e-06 -4.0e-06 11 5682.464 -11342.88 75.32060 4.408400e-17
22 -6.811673 4e-07 -0.0001053 0 -0.0001924 3.00e-06 -0.0000814 NA 3.3e-06 0 NA -4.2e-06 11 5682.288 -11342.53 75.67391 3.694542e-17
14 -6.811598 4e-07 -0.0001083 0 -0.0001751 -1.85e-05 -0.0000989 NA 4.1e-06 0 -3.8e-06 NA 11 5682.250 -11342.46 75.74863 3.559064e-17
21 -6.812095 5e-07 -0.0000908 0 -0.0001927 1.70e-06 0.0000164 0e+00 NA 0 NA -4.1e-06 11 5677.855 -11333.66 84.53872 4.391284e-19
13 -6.812108 5e-07 -0.0000906 0 -0.0001805 -1.91e-05 0.0000197 0e+00 NA 0 -2.7e-06 NA 11 5677.826 -11333.61 84.59689 4.265403e-19
7 -6.811529 6e-07 -0.0001134 0 -0.0001893 -1.57e-05 -0.0001358 -1e-07 6.0e-06 0 NA NA 11 5677.734 -11333.42 84.78076 3.890758e-19
27 -6.811910 5e-07 -0.0001033 0 -0.0001816 -8.20e-06 -0.0000577 0e+00 3.9e-06 NA -2.1e-06 -2.5e-06 12 5676.703 -11329.35 88.85249 5.080069e-20
30 -6.811640 4e-07 -0.0001068 0 -0.0001795 -2.50e-06 -0.0000895 NA 3.7e-06 0 -2.9e-06 -3.4e-06 12 5673.352 -11322.65 95.55317 1.781703e-21
29 -6.812084 4e-07 -0.0000911 0 -0.0001840 -1.30e-06 0.0000157 0e+00 NA 0 -2.0e-06 -3.8e-06 12 5668.947 -11313.84 104.36391 2.175735e-23
23 -6.811572 6e-07 -0.0001115 0 -0.0001901 -4.70e-06 -0.0001226 0e+00 5.3e-06 0 NA -2.4e-06 12 5668.880 -11313.71 104.49749 2.035168e-23
15 -6.811524 6e-07 -0.0001135 0 -0.0001792 -1.72e-05 -0.0001345 0e+00 5.9e-06 0 -2.4e-06 NA 12 5668.812 -11313.57 104.63440 1.900508e-23
31 -6.811560 5e-07 -0.0001119 0 -0.0001812 -7.70e-06 -0.0001236 0e+00 5.4e-06 0 -2.1e-06 -2.0e-06 13 5659.972 -11293.88 124.32315 1.008124e-27

Given that we are running models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we can check to see whether a weighted model isw favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(log(mean.lambda) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(log(mean.mu) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(log(mean.lambda) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(log(mean.mu) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S6: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change the main conclusions drawn in that they did not generate any statistically significant relationships.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Lambda.top gls(model = log(mean.lambda) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~I(1/log(var.lambda)), method = “REML”) 1 8 -16691 -16637 8353
MCC.Lambda.top.unweighted gls(model = log(mean.lambda) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -10674 -10620 5345
anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Mu.top gls(model = log(mean.mu) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~I(1/log(var.mu)), method = “REML”) 1 8 -11418 -11365 5717
MCC.Mu.top.unweighted gls(model = log(mean.mu) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -4024 -3971 2020

We can run the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")
BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")
MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients, confint(MCC.Lambda.top)) %>% tibble::rownames_to_column()

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

BAMM.lambda.plot <-BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients, confint(MCC.Mu.top)) %>% tibble::rownames_to_column()

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

BAMM.mu.plot <-BAMM.mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

PGLS Model Summary

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
)

Figure S8: This figure is the same basic figure as seen in the manuscrip (FIgure 1), it provides model estimates for four response variables across 100 random trees alongside the MCC tree.

From the estimates of 100 trees we can calculate HPD intervals. Note that these do no account for the CI of each estimate. Thus it is a more a representation of tree uncertainty than model uncertainty.

Table S7: Highest posterior density intervals are calculated for the above figure (ie; models using RGB values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -).

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")
DR HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00193 -0.00844 -2.16e-05 -0.00323 -0.00389 -1.4e-06
Upper 0.00156 -0.00182 5.6e-05 0.00479 0.00523 1.53e-06
hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")
ND HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -7.68e-05 -0.00019 -9.72e-07 -9.09e-05 -0.00015 -4.27e-08
Upper 5.97e-05 9.11e-06 1.62e-06 0.000175 0.000117 4.28e-08
hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")
BAMM Speciation HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00637 -0.00453 -0.000144 -0.00606 -0.0109 -3.18e-06
Upper 0.0037 0.00856 5.74e-05 0.0136 0.00699 2.99e-06
hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")
BAMM Extinction HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00539 -0.00747 -0.000134 -0.0137 -0.0118 -3.05e-06
Upper 0.00531 0.0124 0.000122 0.0155 0.0114 3.18e-06

Subsetted analysis with Armenta data

Using the dataset from Armenta et al. (2008) we can conduct a subsetted analysis where dichromatism values are measured from spectrophotmetry and synthesised into a colour discriminability measure, which is a difference measure between the sexes. To make this dataset comprable with ours we use the absolute difference between the sexes. Therefore making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")
MCC.Armenta.model.data$abs.Colour.discriminability <- abs(MCC.Armenta.model.data$Colour.discriminability)

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = abs(Colour.discriminability)))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S9: There is a correlation between the RGB measures and the Spectophotometry measures, The RGB measures seem to be more noisy around the lower values, where s

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(log(mean.lambda) ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(log(mean.mu) ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S8: The dredged models for Models using a subset dataset where sexual dichromatism was measured using spectrophotometry, all four show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Spectrophotometry Dichromatism Dredge Table
(Intercept) abs.Colour.discriminability bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 abs.Colour.discriminability:bioclim4 abs.Colour.discriminability:log(range.size.m2) abs.Colour.discriminability:NPP abs.Colour.discriminability:PC1.LIG abs.Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 0.2077812 -0.0002188 7.0e-06 -0.0000307 5e-07 -0.0008917 0.0010513 NA NA NA NA NA 8 1033.5047 -2050.758 0.00000 9.974822e-01
2 0.2021525 0.0082305 6.6e-06 0.0001692 5e-07 -0.0008672 0.0010498 NA -0.0002958 NA NA NA 9 1027.5845 -2036.854 13.90384 9.543858e-04
8 0.2080321 -0.0000682 6.7e-06 -0.0000426 5e-07 -0.0005900 0.0010626 NA NA NA -0.0004021 NA 9 1027.4340 -2036.553 14.20484 8.210378e-04
16 0.2072009 -0.0003166 7.0e-06 -0.0000055 5e-07 -0.0009221 0.0007620 NA NA NA NA 0.0003857 9 1027.3188 -2036.322 14.43529 7.316795e-04
1 0.2085571 -0.0024081 4.3e-06 -0.0000037 5e-07 -0.0009627 0.0010369 4.80e-06 NA NA NA NA 9 1022.7301 -2027.145 23.61265 7.438452e-06
24 0.2071343 -0.0000863 6.6e-06 -0.0000048 5e-07 -0.0002820 0.0004750 NA NA NA -0.0008968 0.0008019 10 1021.7263 -2023.067 27.69092 9.680505e-07
10 0.2064178 0.0023155 6.7e-06 0.0000149 5e-07 -0.0006057 0.0010613 NA -0.0000838 NA -0.0003719 NA 10 1021.5903 -2022.795 27.96290 8.449648e-07
18 0.2011316 0.0087725 6.7e-06 0.0002101 5e-07 -0.0008964 0.0007539 NA -0.0003183 NA NA 0.0003943 10 1021.4064 -2022.427 28.33081 7.029863e-07
4 0.2090202 -0.0039708 6.6e-06 0.0000038 2e-07 -0.0008146 0.0010507 NA NA 4e-07 NA NA 9 1020.3173 -2022.319 28.43836 6.661846e-07
9 0.2119186 -0.0072833 -3.8e-06 0.0000024 4e-07 0.0004728 0.0010599 1.73e-05 NA NA -0.0021588 NA 10 1018.4257 -2016.465 34.29215 3.568279e-08
3 0.1906298 0.0238655 1.6e-06 0.0006668 4e-07 -0.0009267 0.0010228 7.80e-06 -0.0009674 NA NA NA 10 1017.2485 -2014.111 36.64656 1.099525e-08
17 0.2082477 -0.0021612 4.6e-06 0.0000019 5e-07 -0.0009645 0.0009343 4.20e-06 NA NA NA 0.0001392 10 1016.5624 -2012.739 38.01879 5.536429e-09
26 0.2115099 -0.0066228 6.8e-06 -0.0001604 5e-07 -0.0002217 0.0004455 NA 0.0002299 NA -0.0010072 0.0008469 11 1015.9362 -2009.408 41.34925 1.047198e-09
6 0.2116837 -0.0081513 6.7e-06 -0.0000832 2e-07 -0.0008182 0.0010514 NA 0.0001337 5e-07 NA NA 10 1014.5280 -2008.670 42.08750 7.239709e-10
12 0.2090246 -0.0034064 6.5e-06 -0.0000081 3e-07 -0.0006328 0.0010580 NA NA 4e-07 -0.0002551 NA 10 1014.2134 -2008.041 42.71676 5.285413e-10
20 0.2084568 -0.0043131 6.6e-06 0.0000338 2e-07 -0.0008433 0.0007312 NA NA 4e-07 NA 0.0004259 10 1014.1659 -2007.946 42.81174 5.040264e-10
11 0.1976270 0.0136600 -5.6e-06 0.0005318 4e-07 0.0004453 0.0010479 1.92e-05 -0.0007642 NA -0.0020748 NA 11 1012.8242 -2003.184 47.57326 4.661264e-11
25 0.2110580 -0.0068143 -3.2e-06 0.0000248 4e-07 0.0006083 0.0006664 1.61e-05 NA NA -0.0023731 0.0005373 11 1012.4755 -2002.487 48.27063 3.289058e-11
19 0.1902052 0.0245924 1.3e-06 0.0006857 4e-07 -0.0009247 0.0010743 8.20e-06 -0.0009988 NA NA -0.0000705 11 1011.1192 -1999.774 50.98317 8.473279e-12
5 0.2112711 -0.0105178 1.7e-06 0.0000755 1e-07 -0.0008747 0.0010260 8.10e-06 NA 7e-07 NA NA 10 1010.0432 -1999.700 51.05717 8.165496e-12
14 0.2152129 -0.0129387 6.7e-06 -0.0002139 2e-07 -0.0005837 0.0010617 NA 0.0003106 4e-07 -0.0003358 NA 11 1008.5207 -1994.577 56.18019 6.302792e-13
28 0.2080277 -0.0029321 6.4e-06 0.0000227 3e-07 -0.0003346 0.0005018 NA NA 3e-07 -0.0007456 0.0007600 11 1008.4601 -1994.456 56.30141 5.932118e-13
22 0.2111255 -0.0085019 6.7e-06 -0.0000534 2e-07 -0.0008469 0.0007319 NA 0.0001340 5e-07 NA 0.0004259 11 1008.3773 -1994.291 56.46701 5.460721e-13
13 0.2150226 -0.0163349 -6.9e-06 0.0000891 0e+00 0.0006307 0.0010490 2.14e-05 NA 8e-07 -0.0022521 NA 11 1005.8636 -1989.263 61.49446 4.421320e-14
27 0.2006302 0.0086723 -4.7e-06 0.0004135 4e-07 0.0005527 0.0007595 1.78e-05 -0.0005695 NA -0.0022549 0.0003979 12 1006.8353 -1989.121 61.63628 4.118677e-14
7 0.2018300 0.0039116 8.0e-07 0.0004027 1e-07 -0.0008692 0.0010204 9.10e-06 -0.0004885 6e-07 NA NA 11 1004.3954 -1986.327 64.43070 1.018491e-14
21 0.2113581 -0.0106230 1.6e-06 0.0000747 1e-07 -0.0008738 0.0010493 8.20e-06 NA 7e-07 NA -0.0000318 11 1003.8780 -1985.292 65.46550 6.070883e-15
30 0.2208911 -0.0228961 6.8e-06 -0.0004056 2e-07 -0.0001895 0.0004313 NA 0.0006526 5e-07 -0.0009842 0.0008670 12 1002.8892 -1981.229 69.52846 7.961453e-16
15 0.2117492 -0.0113312 -7.2e-06 0.0002010 0e+00 0.0006174 0.0010469 2.16e-05 -0.0001674 8e-07 -0.0022294 NA 12 1000.1626 -1975.776 74.98178 5.209781e-17
29 0.2142549 -0.0154870 -6.3e-06 0.0000992 0e+00 0.0007138 0.0007815 2.04e-05 NA 8e-07 -0.0023927 0.0003661 12 999.8185 -1975.088 75.66996 3.693025e-17
23 0.2012141 0.0049585 3.0e-07 0.0004324 1e-07 -0.0008651 0.0011085 9.80e-06 -0.0005374 6e-07 NA -0.0001206 12 998.2747 -1972.000 78.75754 7.887190e-18
31 0.2143554 -0.0156416 -6.3e-06 0.0000957 0e+00 0.0007145 0.0007808 2.03e-05 0.0000053 8e-07 -0.0023938 0.0003672 13 994.1598 -1961.678 89.08005 4.522906e-20
kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Spectrophotometry Dichromatism Dredge Table
(Intercept) abs.Colour.discriminability bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 abs.Colour.discriminability:bioclim4 abs.Colour.discriminability:log(range.size.m2) abs.Colour.discriminability:NPP abs.Colour.discriminability:PC1.LIG abs.Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -2.686495 0.0070877 0.0001787 0.0092136 7.8e-06 -0.0316161 0.0155979 NA NA NA NA NA 9 -512.1586 1042.632 0.000000 9.381493e-01
2 -3.007955 0.4954768 0.0001616 0.0207305 7.4e-06 -0.0303962 0.0155527 NA -0.0171356 NA NA NA 10 -514.8054 1049.997 7.364384 2.361111e-02
8 -2.671937 0.0095478 0.0001664 0.0088630 7.5e-06 -0.0204339 0.0165739 NA NA NA -0.0141963 NA 10 -514.9950 1050.376 7.743674 1.953235e-02
16 -2.693478 0.0026019 0.0001826 0.0096006 7.7e-06 -0.0329197 0.0067574 NA NA NA NA 0.0119257 10 -515.1631 1050.712 8.079823 1.651052e-02
24 -2.671555 0.0021899 0.0001618 0.0093363 6.9e-06 -0.0108082 -0.0024601 NA NA NA -0.0302019 0.0271157 11 -516.7325 1055.929 13.296582 1.216062e-03
18 -3.030254 0.5136141 0.0001649 0.0216731 7.3e-06 -0.0317065 0.0062567 NA -0.0179362 NA NA 0.0125392 11 -517.7659 1057.996 15.363318 4.326830e-04
10 -2.888922 0.3319758 0.0001586 0.0165913 7.3e-06 -0.0229029 0.0162654 NA -0.0113381 NA -0.0100322 NA 11 -517.8235 1058.111 15.478518 4.084646e-04
1 -2.684460 -0.0005100 0.0001701 0.0092809 7.8e-06 -0.0319468 0.0154800 0.0000173 NA NA NA NA 10 -520.3848 1061.156 18.523187 8.912774e-05
26 -2.708387 0.0570878 0.0001604 0.0106474 6.9e-06 -0.0113655 -0.0022010 NA -0.0019260 NA -0.0292503 0.0266915 12 -519.6806 1063.911 21.278224 2.247835e-05
4 -2.664734 -0.0766697 0.0001737 0.0101956 2.5e-06 -0.0302954 0.0151253 NA NA 9.40e-06 NA NA 10 -522.3115 1065.009 22.376701 1.297876e-05
9 -2.620545 -0.0880714 0.0000330 0.0091232 7.6e-06 -0.0063902 0.0165039 0.0002312 NA NA -0.0377099 NA 11 -521.8047 1066.073 23.440942 7.623179e-06
3 -3.198180 0.7554332 0.0000956 0.0284381 7.4e-06 -0.0317911 0.0147218 0.0001126 -0.0279927 NA NA NA 11 -522.4668 1067.398 24.765264 3.931546e-06
17 -2.698378 0.0162276 0.0001990 0.0095441 7.6e-06 -0.0324836 0.0056685 -0.0000325 NA NA NA 0.0137191 11 -523.2506 1068.965 26.332758 1.795504e-06
25 -2.624345 -0.0878830 0.0000385 0.0095694 7.1e-06 0.0017317 -0.0015856 0.0002142 NA NA -0.0511968 0.0257625 12 -523.6822 1071.914 29.281408 4.110503e-07
11 -3.103273 0.6223668 -0.0000320 0.0270494 7.3e-06 -0.0071098 0.0157656 0.0003129 -0.0261909 NA -0.0364040 NA 12 -523.9878 1072.525 29.892508 3.028283e-07
6 -2.872767 0.2560291 0.0001646 0.0170736 4.0e-06 -0.0299649 0.0152592 NA -0.0107060 6.30e-06 NA NA 11 -525.1190 1072.702 30.069689 2.771545e-07
12 -2.656835 -0.0592550 0.0001645 0.0097310 3.2e-06 -0.0213688 0.0160145 NA NA 7.70e-06 -0.0116416 NA 11 -525.3025 1073.069 30.436578 2.307029e-07
20 -2.672292 -0.0770310 0.0001777 0.0105126 2.6e-06 -0.0315844 0.0067948 NA NA 9.00e-06 NA 0.0112662 11 -525.3603 1073.185 30.552288 2.177344e-07
19 -3.136655 0.6624291 0.0001249 0.0259779 7.3e-06 -0.0321603 0.0083676 0.0000666 -0.0241379 NA NA 0.0090409 12 -525.5520 1075.653 33.021033 6.336447e-08
28 -2.663124 -0.0357086 0.0001610 0.0097929 4.6e-06 -0.0118446 -0.0017365 NA NA 4.30e-06 -0.0279130 0.0256460 12 -527.2371 1079.024 36.391183 1.174970e-08
27 -2.923316 0.3530486 -0.0000027 0.0206285 6.9e-06 0.0000782 0.0006450 0.0002671 -0.0162510 NA -0.0483668 0.0219478 13 -526.2886 1079.219 36.586890 1.065442e-08
22 -2.914834 0.3101861 0.0001673 0.0185352 4.4e-06 -0.0312776 0.0064477 NA -0.0124597 5.40e-06 NA 0.0119474 12 -528.1204 1080.790 38.157781 4.857527e-09
14 -2.743795 0.0765554 0.0001616 0.0126339 3.8e-06 -0.0221897 0.0159845 NA -0.0044300 6.60e-06 -0.0103857 NA 12 -528.1225 1080.794 38.161874 4.847595e-09
5 -2.646902 -0.1435694 0.0001311 0.0108597 8.0e-07 -0.0314027 0.0143755 0.0000830 NA 1.28e-05 NA NA 11 -530.1941 1082.852 40.219766 1.732450e-09
13 -2.555957 -0.3017598 -0.0000487 0.0112730 -2.1e-06 -0.0003634 0.0152499 0.0003642 NA 1.76e-05 -0.0453855 NA 12 -531.0085 1086.566 43.933901 2.704871e-10
30 -2.603416 -0.1287498 0.0001626 0.0078108 4.2e-06 -0.0111079 -0.0019637 NA 0.0030497 5.00e-06 -0.0290449 0.0260425 13 -530.0513 1086.745 44.112144 2.474239e-10
7 -3.045354 0.4701151 0.0000871 0.0247962 2.7e-06 -0.0314447 0.0141569 0.0001342 -0.0210759 8.90e-06 NA NA 12 -532.6233 1089.796 47.163630 5.380466e-11
21 -2.662261 -0.1088665 0.0001564 0.0107721 1.8e-06 -0.0318399 0.0081816 0.0000396 NA 1.07e-05 NA 0.0089216 12 -533.2753 1091.100 48.467514 2.803399e-11
29 -2.573270 -0.2544878 -0.0000263 0.0111512 -4.0e-07 0.0047402 0.0011884 0.0003213 NA 1.37e-05 -0.0543662 0.0204096 13 -533.4178 1093.478 50.845155 8.538604e-12
15 -2.831379 0.1234432 -0.0000707 0.0207451 -7.0e-07 -0.0017671 0.0150646 0.0003865 -0.0143568 1.47e-05 -0.0433697 NA 13 -533.6497 1093.941 51.309030 6.771077e-12
23 -3.021846 0.4420754 0.0001096 0.0235154 3.3e-06 -0.0317586 0.0096218 0.0000977 -0.0192309 7.60e-06 NA 0.0065712 13 -535.7653 1098.173 55.540280 8.163068e-13
31 -2.716186 -0.0349812 -0.0000391 0.0161109 3.0e-07 0.0036812 0.0019748 0.0003353 -0.0075020 1.24e-05 -0.0527335 0.0191444 14 -536.1666 1101.075 58.442891 1.912315e-13
kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Spectrophotometry Dichromatism Dredge Table
(Intercept) abs.Colour.discriminability bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 abs.Colour.discriminability:bioclim4 abs.Colour.discriminability:log(range.size.m2) abs.Colour.discriminability:NPP abs.Colour.discriminability:PC1.LIG abs.Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -3.102033 0.0005342 -3.25e-05 0.0016996 -1.4e-06 -0.0019005 0.0021590 NA NA NA NA NA 8 254.4471 -492.6424 0.000000 9.865856e-01
2 -3.011015 -0.1423896 -2.74e-05 -0.0016071 -1.3e-06 -0.0023624 0.0022014 NA 0.0049897 NA NA NA 9 250.6065 -482.8977 9.744703 7.552653e-03
8 -3.102034 0.0005328 -3.25e-05 0.0016996 -1.4e-06 -0.0019034 0.0021588 NA NA NA 0.0000038 NA 9 249.6525 -480.9898 11.652618 2.909383e-03
16 -3.100824 0.0011510 -3.30e-05 0.0016302 -1.4e-06 -0.0017539 0.0032438 NA NA NA NA -0.0014356 9 249.6368 -480.9583 11.684104 2.863940e-03
10 -2.984592 -0.1824367 -2.76e-05 -0.0025773 -1.3e-06 -0.0006878 0.0023634 NA 0.0064189 NA -0.0024181 NA 10 246.1462 -471.9064 20.735999 3.100059e-05
18 -3.007807 -0.1446367 -2.77e-05 -0.0017519 -1.3e-06 -0.0022091 0.0034060 NA 0.0050920 NA NA -0.0015931 10 245.8256 -471.2652 21.377207 2.249745e-05
1 -3.102052 0.0006060 -3.25e-05 0.0016988 -1.4e-06 -0.0018980 0.0021598 -2.00e-07 NA NA NA NA 9 244.7581 -471.2009 21.441484 2.178591e-05
24 -3.100506 0.0009393 -3.22e-05 0.0016134 -1.3e-06 -0.0026428 0.0036166 NA NA NA 0.0012712 -0.0020332 10 245.0686 -469.7511 22.891224 1.055282e-05
4 -3.101388 -0.0013224 -3.27e-05 0.0017160 -1.5e-06 -0.0018653 0.0021526 NA NA 2.0e-07 NA NA 9 242.3609 -466.4066 26.235779 1.982019e-06
3 -2.975150 -0.1922942 -1.26e-05 -0.0031422 -1.3e-06 -0.0021749 0.0023320 -2.33e-05 0.0071143 NA NA NA 10 241.4066 -462.4272 30.215135 2.710202e-07
26 -2.989499 -0.1737415 -2.78e-05 -0.0024101 -1.3e-06 -0.0010229 0.0029027 NA 0.0061192 NA -0.0018279 -0.0007655 11 241.4495 -460.4351 32.207318 1.000932e-07
9 -3.102090 0.0007097 -3.23e-05 0.0016978 -1.4e-06 -0.0019276 0.0021582 -4.00e-07 NA NA 0.0000452 NA 10 240.3903 -460.3947 32.247714 9.809181e-08
17 -3.099659 -0.0019395 -3.68e-05 0.0016502 -1.4e-06 -0.0018305 0.0034866 6.90e-06 NA NA NA -0.0018010 10 240.0956 -459.8052 32.837151 7.305301e-08
6 -2.965005 -0.2236032 -2.65e-05 -0.0028552 -2.4e-06 -0.0021865 0.0021511 NA 0.0071363 2.2e-06 NA NA 10 239.0351 -457.6843 34.958039 2.529840e-08
12 -3.101362 -0.0014584 -3.26e-05 0.0017177 -1.5e-06 -0.0019218 0.0021473 NA NA 2.0e-07 0.0000783 NA 10 237.5941 -454.8022 37.840122 5.987653e-09
20 -3.100049 -0.0010492 -3.31e-05 0.0016492 -1.5e-06 -0.0017109 0.0032446 NA NA 2.0e-07 NA -0.0014467 10 237.5537 -454.7213 37.921034 5.750252e-09
11 -2.973061 -0.1962325 -1.53e-05 -0.0031810 -1.3e-06 -0.0016808 0.0023603 -1.93e-05 0.0071952 NA -0.0007606 NA 11 237.0550 -451.6461 40.996266 1.235687e-09
19 -2.978672 -0.1865916 -1.47e-05 -0.0029942 -1.3e-06 -0.0021440 0.0027509 -2.03e-05 0.0068739 NA NA -0.0005765 11 236.6406 -450.8172 41.825217 8.164026e-10
25 -3.100294 0.0003062 -3.31e-05 0.0016195 -1.3e-06 -0.0025599 0.0036266 1.50e-06 NA NA 0.0011294 -0.0020442 11 235.8111 -449.1582 43.484205 3.561719e-10
14 -2.936793 -0.2665863 -2.68e-05 -0.0038823 -2.5e-06 -0.0004524 0.0023176 NA 0.0086555 2.3e-06 -0.0024991 NA 11 234.5894 -446.7148 45.927584 1.049752e-10
22 -2.959170 -0.2301133 -2.69e-05 -0.0030776 -2.5e-06 -0.0020078 0.0034852 NA 0.0073560 2.3e-06 NA -0.0017674 11 234.2914 -446.1189 46.523451 7.792839e-11
5 -3.101147 -0.0021301 -3.32e-05 0.0017245 -1.5e-06 -0.0018753 0.0021466 1.00e-06 NA 2.0e-07 NA NA 10 232.7712 -445.1564 47.486009 4.815916e-11
28 -3.098868 -0.0035688 -3.24e-05 0.0016485 -1.6e-06 -0.0027373 0.0036944 NA NA 5.0e-07 0.0015312 -0.0021780 11 233.0483 -443.6326 49.009797 2.247980e-11
27 -2.976536 -0.1904124 -1.60e-05 -0.0030519 -1.3e-06 -0.0018275 0.0026779 -1.83e-05 0.0069793 NA -0.0004974 -0.0004506 12 232.3627 -440.1761 52.466258 3.992378e-12
7 -2.942202 -0.2527996 -1.45e-05 -0.0039375 -2.2e-06 -0.0020578 0.0022663 -1.92e-05 0.0085725 1.9e-06 NA NA 11 229.7460 -437.0281 55.614256 8.272802e-13
30 -2.942035 -0.2572896 -2.69e-05 -0.0036970 -2.5e-06 -0.0008737 0.0030026 NA 0.0083212 2.3e-06 -0.0017500 -0.0009736 12 229.9093 -435.2692 57.373161 3.433295e-13
13 -3.101109 -0.0022357 -3.34e-05 0.0017254 -1.5e-06 -0.0018530 0.0021476 1.20e-06 NA 3.0e-07 -0.0000338 NA 11 228.4163 -434.3686 58.273768 2.188501e-13
21 -3.096804 -0.0099118 -3.95e-05 0.0017149 -1.7e-06 -0.0017578 0.0036361 1.11e-05 NA 7.0e-07 NA -0.0020538 11 228.1812 -433.8984 58.743975 1.729986e-13
15 -2.934658 -0.2668675 -2.02e-05 -0.0040951 -2.4e-06 -0.0010031 0.0023198 -1.03e-05 0.0088820 2.1e-06 -0.0016061 NA 12 225.4547 -426.3600 66.282335 3.991248e-15
23 -2.945869 -0.2475112 -1.85e-05 -0.0037298 -2.3e-06 -0.0019886 0.0030527 -1.31e-05 0.0082503 2.1e-06 NA -0.0010905 12 225.0366 -425.5238 67.118540 2.627417e-15
29 -3.097523 -0.0074677 -3.60e-05 0.0016836 -1.7e-06 -0.0024325 0.0037579 5.90e-06 NA 7.0e-07 0.0010393 -0.0022656 12 223.8961 -423.2430 69.399375 8.399473e-16
31 -2.939504 -0.2589943 -2.17e-05 -0.0038954 -2.4e-06 -0.0012438 0.0028952 -8.20e-06 0.0085543 2.2e-06 -0.0011604 -0.0008185 13 220.7883 -414.9347 77.707705 1.318624e-17
kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Spectrophotometry Dichromatism Dredge Table
(Intercept) abs.Colour.discriminability bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 abs.Colour.discriminability:bioclim4 abs.Colour.discriminability:log(range.size.m2) abs.Colour.discriminability:NPP abs.Colour.discriminability:PC1.LIG abs.Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -6.792943 -0.0103858 3.41e-05 -0.0002965 -8.0e-07 -0.0106156 0.0000352 NA NA NA NA NA 8 -80.41697 177.0857 0.000000 9.755486e-01
2 -6.629991 -0.2696903 4.33e-05 -0.0062547 -6.0e-07 -0.0114938 0.0000244 NA 0.0090568 NA NA NA 9 -83.65433 185.6239 8.538207 1.365214e-02
16 -6.791609 -0.0089665 3.27e-05 -0.0004075 -8.0e-07 -0.0102634 0.0025209 NA NA NA NA -0.0033232 9 -84.56739 187.4500 10.364333 5.478512e-03
8 -6.792929 -0.0103492 3.41e-05 -0.0002974 -8.0e-07 -0.0105457 0.0000407 NA NA NA -0.0000938 NA 9 -84.64806 187.6114 10.525683 5.053890e-03
10 -6.578647 -0.3485384 4.31e-05 -0.0081510 -6.0e-07 -0.0083501 0.0002912 NA 0.0118731 NA -0.0045846 NA 10 -87.51374 195.4134 18.327755 1.021945e-04
18 -6.619830 -0.2817268 4.22e-05 -0.0066989 -6.0e-07 -0.0111292 0.0029242 NA 0.0095350 NA NA -0.0038774 10 -87.72589 195.8377 18.752067 8.265882e-05
1 -6.793837 -0.0071012 3.79e-05 -0.0003277 -8.0e-07 -0.0104952 0.0000638 -7.00e-06 NA NA NA NA 9 -89.52020 197.3556 20.269961 3.869749e-05
24 -6.791457 -0.0095030 3.41e-05 -0.0004282 -7.0e-07 -0.0123313 0.0034424 NA NA NA 0.0029834 -0.0047902 10 -88.51791 197.4218 20.336105 3.743862e-05
4 -6.792362 -0.0119234 3.40e-05 -0.0002822 -9.0e-07 -0.0105858 0.0000320 NA NA 2.0e-07 NA NA 9 -91.92854 202.1723 25.086646 3.481392e-06
3 -6.546620 -0.3881417 7.70e-05 -0.0097911 -5.0e-07 -0.0110606 0.0002367 -5.34e-05 0.0140699 NA NA NA 10 -92.07196 204.5299 27.444213 1.071060e-06
26 -6.595359 -0.3203397 4.24e-05 -0.0075827 -6.0e-07 -0.0095563 0.0022257 NA 0.0109002 NA -0.0024481 -0.0027530 11 -91.56060 205.5852 28.499491 6.319213e-07
9 -6.795188 -0.0035460 4.36e-05 -0.0003527 -8.0e-07 -0.0114750 0.0000104 -1.58e-05 NA NA 0.0015187 NA 10 -93.28777 206.9615 29.875824 3.175383e-07
17 -6.790340 -0.0127444 2.80e-05 -0.0003852 -8.0e-07 -0.0103599 0.0028331 8.50e-06 NA NA NA -0.0037868 10 -93.54185 207.4697 30.383992 2.462912e-07
6 -6.545265 -0.4185973 4.46e-05 -0.0085371 -2.7e-06 -0.0111960 -0.0000539 NA 0.0130245 4.0e-06 NA NA 10 -94.64428 209.6745 32.588849 8.178437e-08
20 -6.791009 -0.0105540 3.26e-05 -0.0003928 -9.0e-07 -0.0102326 0.0025179 NA NA 2.0e-07 NA -0.0033237 10 -96.07841 212.5428 35.457096 1.949119e-08
12 -6.792376 -0.0118581 3.40e-05 -0.0002830 -9.0e-07 -0.0105609 0.0000341 NA NA 2.0e-07 -0.0000348 NA 10 -96.12841 212.6428 35.557100 1.854055e-08
11 -6.546790 -0.3878181 7.72e-05 -0.0097879 -5.0e-07 -0.0111019 0.0002344 -5.38e-05 0.0140635 NA 0.0000644 NA 11 -95.84961 214.1632 37.077509 8.669023e-09
19 -6.556102 -0.3736380 7.09e-05 -0.0093945 -5.0e-07 -0.0109759 0.0014473 -4.46e-05 0.0134426 NA NA -0.0016659 11 -96.24076 214.9455 37.859814 5.862662e-09
25 -6.793167 -0.0044397 4.12e-05 -0.0004673 -7.0e-07 -0.0129949 0.0033627 -1.18e-05 NA NA 0.0041352 -0.0047097 11 -97.16621 216.7964 39.710712 2.323685e-09
14 -6.490797 -0.5026104 4.45e-05 -0.0105322 -2.7e-06 -0.0079572 0.0002189 NA 0.0159988 4.0e-06 -0.0047149 NA 11 -98.48958 219.4431 42.357454 6.186494e-10
22 -6.529743 -0.4396690 4.35e-05 -0.0091370 -2.8e-06 -0.0107908 0.0030334 NA 0.0137539 4.2e-06 NA -0.0041341 11 -98.67549 219.8150 42.729279 5.136933e-10
5 -6.794334 -0.0057058 3.83e-05 -0.0003407 -7.0e-07 -0.0105065 0.0000685 -7.60e-06 NA -1.0e-07 NA NA 10 -100.92841 222.2428 45.157109 1.525834e-10
28 -6.788932 -0.0162012 3.38e-05 -0.0003694 -1.1e-06 -0.0124711 0.0035496 NA NA 7.0e-07 0.0033702 -0.0049825 11 -99.96767 222.3993 45.313628 1.410976e-10
27 -6.561226 -0.3645445 7.41e-05 -0.0092555 -5.0e-07 -0.0117654 0.0016335 -4.98e-05 0.0131979 NA 0.0012566 -0.0019849 12 -99.92687 224.4030 47.317360 5.181015e-11
7 -6.489092 -0.4927366 7.38e-05 -0.0111805 -2.2e-06 -0.0108759 0.0001476 -4.68e-05 0.0166261 3.2e-06 NA NA 11 -103.17758 228.8191 51.733441 5.694752e-12
30 -6.506666 -0.4760139 4.37e-05 -0.0099741 -2.8e-06 -0.0092833 0.0023628 NA 0.0150450 4.2e-06 -0.0023486 -0.0030542 12 -102.51344 229.5762 52.490494 3.900164e-12
13 -6.796232 -0.0006312 4.47e-05 -0.0003795 -6.0e-07 -0.0115476 0.0000170 -1.75e-05 NA -2.0e-07 0.0015968 NA 11 -104.68194 231.8278 54.742163 1.265142e-12
21 -6.787454 -0.0205873 2.53e-05 -0.0003187 -1.1e-06 -0.0102908 0.0029741 1.26e-05 NA 7.0e-07 NA -0.0040147 11 -104.90431 232.2726 55.186921 1.012888e-12
15 -6.483098 -0.5039427 6.92e-05 -0.0113082 -2.3e-06 -0.0100458 0.0001901 -3.95e-05 0.0168697 3.3e-06 -0.0012815 NA 12 -106.92194 238.3932 61.307488 4.747851e-14
23 -6.496504 -0.4833026 6.45e-05 -0.0107540 -2.4e-06 -0.0107296 0.0019254 -3.29e-05 0.0159927 3.5e-06 NA -0.0024605 12 -107.28539 239.1201 62.034394 3.301045e-14
29 -6.790799 -0.0108871 3.88e-05 -0.0004124 -1.0e-06 -0.0128884 0.0034655 -8.10e-06 NA 5.0e-07 0.0040550 -0.0048745 12 -108.52941 241.6081 64.522435 9.514371e-15
31 -6.497205 -0.4820471 6.49e-05 -0.0107357 -2.4e-06 -0.0108121 0.0019433 -3.35e-05 0.0159601 3.5e-06 0.0001301 -0.0024913 13 -110.96642 248.5748 71.489128 2.921340e-16

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

Armenta.data$abs.Colour.discriminability <- abs(Armenta.data$Colour.discriminability)

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(abs.Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(log(mean.lambda) ~ abs.Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(log(mean.mu) ~ abs.Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ abs.Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Supp Fig S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `abs.Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

####################################
#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

######################

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

######################

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

#########################################

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S9: Measures of sexual dichromatism using a more precise measure in spectrophotometry do not change the main patterns seen within the PGLS models analysed here. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error; additionally these environmental predictors were not measured differently in this analysis compared to the previous above so we have no obvious reason to make conclusions from this subset for the effects of the environmental predictors.

Table S9: Highest posterior density intervals are calculated for the above figure (ie; models using Spectrophotometry values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

###############
Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

################
Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")
Armenta DR HPD Interval
  abs.Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0393 -0.0143 -0.000139 0.0108 -0.0427 -4.27e-06
Upper 0.0421 0.0316 0.000159 0.0347 0.0076 1.44e-05
Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")
Armenta ND HPD Interval
  abs.Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0021 -0.000828 -4.96e-06 0.000491 -0.00228 -6.59e-08
Upper 0.0026 0.00121 1.38e-05 0.00181 0.000572 6.88e-07
Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")
Armenta BAMM Speciation HPD Interval
  abs.Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0243 -0.00942 -6.03e-05 -0.00231 -0.0114 -3.9e-06
Upper 0.0144 0.00725 2.65e-05 0.00543 0.00317 3.2e-06
Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")
Armenta BAMM Extinction HPD Interval
  abs.Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0337 -0.0114 -5.71e-05 -0.00565 -0.0228 -5.39e-06
Upper 0.0269 0.0112 7.6e-05 0.00487 0.00156 4.43e-06

Analysis using male-bias measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism, and paternal care) is expected to be relatively noisy given the precision of the measurment used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species within the orginal dataset of sexual dichromatism scores from RGB measures. This male-bias sexual selection score is based on three components, combined in a phylogenetic pca (ppca) with the following loadings:

  • sexual size dimorphism: 0.37
  • social polygyny: 0.57
  • paternal care: -0.57

As such, high values for this score are expected to have high sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset are calculated in the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

Figure S10: The relationship between absolute sexual dichromatism and the sexual selection. Although positively correlated the distribution is very noisy with species having high sexual dichromatism and low sexual selection and visa versa.

#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(log(mean.lambda) ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(log(mean.mu) ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S10: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms are not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 0.1484807 -5e-07 -0.0001800 0 0.0001073 0.0000788 0.0004381 NA NA NA NA NA 8 5451.548 -10887.04 0.00000 9.994587e-01
16 0.1482246 -5e-07 -0.0001720 0 0.0000929 0.0000644 0.0003782 NA NA NA NA -0.0001167 9 5444.064 -10870.05 16.98341 2.050523e-04
8 0.1479484 -5e-07 -0.0001570 0 0.0000777 0.0001046 0.0004928 NA NA NA 0.0001280 NA 9 5444.026 -10869.98 17.05889 1.974574e-04
2 0.1483243 -5e-07 -0.0001765 0 0.0001017 0.0000834 0.0018018 NA -4.87e-05 NA NA NA 9 5443.645 -10869.22 17.82105 1.348878e-04
1 0.1478764 -9e-07 -0.0001561 0 0.0000738 0.0000871 0.0012068 -2.3e-06 NA NA NA NA 9 5440.065 -10862.06 24.98016 3.761776e-06
4 0.1484220 -6e-07 -0.0001730 0 0.0001024 0.0000814 0.0008188 NA NA 0e+00 NA NA 9 5436.296 -10854.52 32.51872 8.677880e-08
24 0.1477981 -5e-07 -0.0001528 0 0.0000692 0.0000925 0.0004458 NA NA NA 0.0001196 -0.0000845 10 5436.455 -10852.82 34.21679 3.712640e-08
18 0.1479767 -6e-07 -0.0001662 0 0.0000832 0.0000688 0.0022340 NA -6.65e-05 NA NA -0.0001323 10 5436.226 -10852.36 34.67399 2.953948e-08
10 0.1478037 -5e-07 -0.0001538 0 0.0000725 0.0001088 0.0017783 NA -4.59e-05 NA 0.0001274 NA 10 5436.116 -10852.14 34.89512 2.644759e-08
17 0.1478092 -9e-07 -0.0001543 0 0.0000700 0.0000806 0.0011337 -2.2e-06 NA NA NA -0.0000480 10 5432.459 -10844.83 42.20969 6.824126e-10
9 0.1474813 -8e-07 -0.0001389 0 0.0000518 0.0001085 0.0011728 -2.1e-06 NA NA 0.0001104 NA 10 5432.373 -10844.66 42.38150 6.262388e-10
3 0.1478909 -9e-07 -0.0001563 0 0.0000742 0.0000866 0.0010291 -2.3e-06 6.70e-06 NA NA NA 10 5432.150 -10844.21 42.82661 5.012840e-10
20 0.1482076 -6e-07 -0.0001673 0 0.0000906 0.0000682 0.0006967 NA NA 0e+00 NA -0.0001026 10 5428.776 -10837.46 49.57473 1.716915e-11
12 0.1478969 -6e-07 -0.0001505 0 0.0000733 0.0001068 0.0008569 NA NA 0e+00 0.0001269 NA 10 5428.762 -10837.43 49.60299 1.692825e-11
6 0.1482266 -6e-07 -0.0001680 0 0.0000951 0.0000872 0.0025093 NA -5.90e-05 0e+00 NA NA 10 5428.428 -10836.77 50.27032 1.212563e-11
26 0.1475846 -6e-07 -0.0001479 0 0.0000610 0.0000958 0.0021036 NA -5.95e-05 NA 0.0001173 -0.0000990 11 5428.596 -10835.09 51.95195 5.230492e-12
5 0.1477933 -1e-06 -0.0001475 0 0.0000674 0.0000903 0.0016611 -2.4e-06 NA -1e-07 NA NA 10 5424.857 -10829.62 57.41215 3.410950e-13
25 0.1474541 -8e-07 -0.0001383 0 0.0000503 0.0001049 0.0011368 -2.0e-06 NA NA 0.0001086 -0.0000240 11 5424.751 -10827.39 59.64290 1.118085e-13
19 0.1477962 -9e-07 -0.0001541 0 0.0000696 0.0000808 0.0012568 -2.1e-06 -4.70e-06 NA NA -0.0000499 11 5424.591 -10827.07 59.96353 9.524668e-14
11 0.1474880 -8e-07 -0.0001390 0 0.0000520 0.0001082 0.0010942 -2.1e-06 2.90e-06 NA 0.0001103 NA 11 5424.457 -10826.81 60.23059 8.334079e-14
28 0.1477791 -6e-07 -0.0001479 0 0.0000668 0.0000965 0.0007717 NA NA 0e+00 0.0001201 -0.0000699 11 5421.172 -10820.24 66.80118 3.119414e-15
22 0.1479328 -7e-07 -0.0001603 0 0.0000797 0.0000735 0.0027753 NA -7.32e-05 0e+00 NA -0.0001181 11 5420.966 -10819.82 67.21289 2.539046e-15
14 0.1477158 -6e-07 -0.0001460 0 0.0000665 0.0001122 0.0024557 NA -5.58e-05 0e+00 0.0001260 NA 11 5420.885 -10819.66 67.37476 2.341650e-15
21 0.1477617 -1e-06 -0.0001470 0 0.0000656 0.0000868 0.0016033 -2.3e-06 NA 0e+00 NA -0.0000251 11 5417.249 -10812.39 74.64696 6.171400e-17
13 0.1474084 -1e-06 -0.0001309 0 0.0000460 0.0001112 0.0016076 -2.1e-06 NA 0e+00 0.0001085 NA 11 5417.149 -10812.19 74.84702 5.583943e-17
7 0.1477851 -1e-06 -0.0001473 0 0.0000671 0.0000906 0.0017584 -2.4e-06 -3.60e-06 -1e-07 NA NA 11 5416.948 -10811.79 75.24841 4.568560e-17
27 0.1474468 -8e-07 -0.0001381 0 0.0000501 0.0001050 0.0012077 -2.0e-06 -2.70e-06 NA 0.0001085 -0.0000251 12 5416.883 -10809.64 77.39880 1.558939e-17
30 0.1475395 -6e-07 -0.0001419 0 0.0000575 0.0001007 0.0026498 NA -6.62e-05 0e+00 0.0001176 -0.0000846 12 5413.338 -10802.55 84.48797 4.502300e-19
29 0.1474072 -1e-06 -0.0001309 0 0.0000459 0.0001110 0.0016049 -2.1e-06 NA 0e+00 0.0001084 -0.0000012 12 5409.540 -10794.95 92.08377 1.009318e-20
23 0.1477348 -1e-06 -0.0001465 0 0.0000647 0.0000871 0.0018596 -2.2e-06 -9.70e-06 0e+00 NA -0.0000288 12 5409.384 -10794.64 92.39589 8.634816e-21
15 0.1473924 -1e-06 -0.0001306 0 0.0000455 0.0001118 0.0017935 -2.1e-06 -6.80e-06 0e+00 0.0001087 NA 12 5409.241 -10794.35 92.68281 7.480788e-21
31 0.1473862 -1e-06 -0.0001305 0 0.0000452 0.0001112 0.0018080 -2.1e-06 -7.70e-06 0e+00 0.0001084 -0.0000042 13 5401.675 -10777.20 109.83562 1.410139e-24
kable(Subset.dredged.DR.model, "html", caption = "DR Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.644369 2.45e-05 -0.0113055 -7.0e-07 0.0068813 0.0004409 0.0388722 NA NA NA NA NA 9 -1717.437 3452.947 0.00000 9.851357e-01
16 -2.639661 2.53e-05 -0.0115625 -7.0e-07 0.0066554 -0.0000035 0.0381221 NA NA NA NA -0.0043664 10 -1721.450 3462.989 10.04197 6.499963e-03
2 -2.651446 2.27e-05 -0.0111253 -7.0e-07 0.0068731 0.0003113 0.1032726 NA -0.0023191 NA NA NA 10 -1721.901 3463.892 10.94493 4.138415e-03
8 -2.643685 2.47e-05 -0.0113672 -8.0e-07 0.0067933 0.0002826 0.0391218 NA NA NA -0.0017643 NA 10 -1721.915 3463.919 10.97215 4.082471e-03
1 -2.644777 1.73e-05 -0.0112474 -8.0e-07 0.0069016 0.0002884 0.0537640 -4.40e-05 NA NA NA NA 10 -1726.048 3472.186 19.23901 6.543319e-05
18 -2.648357 2.30e-05 -0.0113613 -7.0e-07 0.0066118 -0.0002393 0.1227636 NA -0.0030524 NA NA -0.0050116 11 -1725.773 3473.654 20.70712 3.140527e-05
24 -2.639316 2.54e-05 -0.0115972 -7.0e-07 0.0066022 -0.0000984 0.0383393 NA NA NA -0.0012567 -0.0041819 11 -1725.963 3474.034 21.08668 2.597662e-05
10 -2.650479 2.30e-05 -0.0111903 -7.0e-07 0.0067936 0.0001735 0.1006675 NA -0.0022170 NA -0.0015983 NA 11 -1726.392 3474.892 21.94541 1.690876e-05
4 -2.639136 2.37e-05 -0.0113014 -1.1e-06 0.0067364 0.0005491 0.0550570 NA NA -1.6e-06 NA NA 10 -1729.238 3478.565 25.61792 2.695472e-06
17 -2.641237 1.95e-05 -0.0114476 -8.0e-07 0.0067343 0.0000040 0.0497120 -3.36e-05 NA NA NA -0.0031418 11 -1730.243 3482.593 29.64645 3.596263e-07
9 -2.644302 1.77e-05 -0.0112841 -8.0e-07 0.0068566 0.0002153 0.0532073 -4.20e-05 NA NA -0.0008810 NA 11 -1730.569 3483.246 30.29926 2.594752e-07
3 -2.648526 1.70e-05 -0.0111570 -7.0e-07 0.0068950 0.0002315 0.0876844 -3.97e-05 -0.0012746 NA NA NA 11 -1730.598 3483.303 30.35634 2.521738e-07
26 -2.647793 2.32e-05 -0.0113950 -7.0e-07 0.0065727 -0.0003052 0.1206409 NA -0.0029697 NA -0.0009490 -0.0048550 12 -1730.301 3484.728 31.78144 1.236643e-07
20 -2.635570 2.45e-05 -0.0115406 -1.0e-06 0.0065504 0.0001210 0.0518306 NA NA -1.4e-06 NA -0.0040359 11 -1733.315 3488.737 35.78960 1.666795e-08
6 -2.646854 2.14e-05 -0.0110830 -1.1e-06 0.0067038 0.0004113 0.1349625 NA -0.0027849 -1.9e-06 NA NA 11 -1733.615 3489.338 36.39048 1.234251e-08
12 -2.638336 2.39e-05 -0.0113656 -1.1e-06 0.0066414 0.0003852 0.0556198 NA NA -1.7e-06 -0.0018515 NA 11 -1733.706 3489.519 36.57214 1.127085e-08
19 -2.647224 1.97e-05 -0.0113350 -7.0e-07 0.0066756 -0.0001778 0.1100270 -2.27e-05 -0.0023104 NA NA -0.0040281 12 -1734.642 3493.410 40.46326 1.610686e-09
25 -2.640864 1.99e-05 -0.0114754 -8.0e-07 0.0067004 -0.0000514 0.0493311 -3.21e-05 NA NA -0.0007139 -0.0030923 12 -1734.769 3493.665 40.71787 1.418151e-09
11 -2.648042 1.75e-05 -0.0111939 -8.0e-07 0.0068500 0.0001586 0.0870861 -3.76e-05 -0.0012730 NA -0.0008797 NA 12 -1735.118 3494.364 41.41693 9.998225e-10
5 -2.639640 1.65e-05 -0.0112442 -1.1e-06 0.0067596 0.0003957 0.0694708 -4.36e-05 NA -1.6e-06 NA NA 11 -1737.859 3497.825 44.87789 1.771687e-10
22 -2.644516 2.18e-05 -0.0113100 -1.1e-06 0.0064809 -0.0001171 0.1492103 NA -0.0034126 -1.7e-06 NA -0.0046899 12 -1737.559 3499.245 46.29775 8.710986e-11
28 -2.635111 2.47e-05 -0.0115772 -1.1e-06 0.0064896 0.0000208 0.0524230 NA NA -1.4e-06 -0.0013747 -0.0038249 12 -1737.818 3499.763 46.81566 6.723639e-11
14 -2.645817 2.17e-05 -0.0111498 -1.2e-06 0.0066197 0.0002691 0.1325514 NA -0.0026832 -1.9e-06 -0.0016616 NA 12 -1738.100 3500.327 47.37978 5.071161e-11
27 -2.646817 2.00e-05 -0.0113622 -7.0e-07 0.0066447 -0.0002281 0.1093216 -2.14e-05 -0.0022967 NA -0.0006587 -0.0039771 13 -1739.168 3504.485 51.53820 6.340417e-12
21 -2.637030 1.86e-05 -0.0114218 -1.1e-06 0.0066273 0.0001336 0.0642539 -3.45e-05 NA -1.4e-06 NA -0.0027665 12 -1742.094 3508.315 55.36799 9.343053e-13
7 -2.644376 1.61e-05 -0.0111162 -1.1e-06 0.0067357 0.0003290 0.1181574 -3.75e-05 -0.0017672 -1.8e-06 NA NA 12 -1742.351 3508.828 55.88132 7.228072e-13
13 -2.639083 1.70e-05 -0.0112841 -1.1e-06 0.0067082 0.0003160 0.0690252 -4.13e-05 NA -1.6e-06 -0.0009791 NA 12 -1742.374 3508.875 55.92791 7.061633e-13
30 -2.643870 2.20e-05 -0.0113459 -1.1e-06 0.0064361 -0.0001880 0.1471758 NA -0.0033252 -1.7e-06 -0.0010505 -0.0045120 13 -1742.080 3510.308 57.36050 3.450013e-13
23 -2.643470 1.87e-05 -0.0112857 -1.1e-06 0.0065443 -0.0000596 0.1365454 -2.19e-05 -0.0026920 -1.6e-06 NA -0.0037463 13 -1746.434 3519.016 66.06877 4.434500e-15
29 -2.636582 1.90e-05 -0.0114524 -1.1e-06 0.0065870 0.0000719 0.0640021 -3.28e-05 NA -1.4e-06 -0.0008236 -0.0027038 13 -1746.614 3519.376 66.42937 3.702894e-15
15 -2.643820 1.66e-05 -0.0111564 -1.1e-06 0.0066838 0.0002487 0.1177972 -3.52e-05 -0.0017704 -1.8e-06 -0.0009854 NA 13 -1746.865 3519.878 66.93090 2.881620e-15
31 -2.642995 1.90e-05 -0.0113159 -1.1e-06 0.0065069 -0.0001168 0.1359737 -2.04e-05 -0.0026794 -1.6e-06 -0.0007727 -0.0036828 14 -1750.955 3530.081 77.13422 1.753934e-17
kable(Subset.dredged.spec.model, "html", caption = "BAMM Speciation Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.835512 -7.3e-06 0.0004994 0 -0.0003236 0.0004112 -0.0003583 NA NA NA NA NA 8 2617.018 -5217.977 0.00000 9.982022e-01
8 -2.837029 -7.4e-06 0.0005575 0 -0.0003517 0.0004687 -0.0004496 NA NA NA 0.0004539 NA 9 2610.740 -5203.407 14.56967 6.846268e-04
16 -2.835892 -7.3e-06 0.0005110 0 -0.0003416 0.0003737 -0.0004570 NA NA NA NA -0.0002267 9 2610.652 -5203.231 14.74603 6.268432e-04
2 -2.835374 -7.2e-06 0.0004942 0 -0.0003057 0.0003976 -0.0048344 NA 0.0001612 NA NA NA 9 2610.384 -5202.694 15.28266 4.793257e-04
1 -2.835131 -7.0e-06 0.0004827 0 -0.0003044 0.0004030 -0.0008253 1.5e-06 NA NA NA NA 9 2605.952 -5193.830 24.14720 5.697974e-06
24 -2.837371 -7.4e-06 0.0005678 0 -0.0003683 0.0004333 -0.0005410 NA NA NA 0.0004501 -0.0002114 10 2604.367 -5188.644 29.33283 4.262599e-07
10 -2.836889 -7.3e-06 0.0005522 0 -0.0003347 0.0004556 -0.0046804 NA 0.0001523 NA 0.0004513 NA 10 2604.101 -5188.112 29.86480 3.267087e-07
18 -2.835724 -7.3e-06 0.0005051 0 -0.0003243 0.0003679 -0.0041362 NA 0.0001330 NA NA -0.0001943 10 2604.018 -5187.947 30.03002 3.008033e-07
4 -2.835551 -7.3e-06 0.0005019 0 -0.0003257 0.0004112 -0.0001576 NA NA 0 NA NA 9 2602.884 -5187.695 30.28165 2.652418e-07
17 -2.835405 -6.9e-06 0.0004876 0 -0.0003165 0.0003449 -0.0012820 2.5e-06 NA NA NA -0.0003181 10 2599.694 -5179.298 38.67931 3.982116e-09
9 -2.836626 -7.1e-06 0.0005397 0 -0.0003311 0.0004602 -0.0009541 1.6e-06 NA NA 0.0004565 NA 10 2599.679 -5179.269 38.70807 3.925262e-09
3 -2.835149 -7.1e-06 0.0004842 0 -0.0002961 0.0003945 -0.0044284 1.0e-06 0.0001357 NA NA NA 10 2599.353 -5178.616 39.36092 2.832071e-09
26 -2.837205 -7.4e-06 0.0005620 0 -0.0003517 0.0004276 -0.0040318 NA 0.0001262 NA 0.0004485 -0.0001808 11 2597.730 -5173.352 44.62511 2.037042e-10
12 -2.837067 -7.4e-06 0.0005600 0 -0.0003538 0.0004688 -0.0002498 NA NA 0 0.0004539 NA 10 2596.607 -5173.124 44.85298 1.817677e-10
20 -2.835901 -7.3e-06 0.0005117 0 -0.0003421 0.0003741 -0.0003955 NA NA 0 NA -0.0002247 10 2596.527 -5172.964 45.01281 1.678074e-10
6 -2.835390 -7.3e-06 0.0004952 0 -0.0003067 0.0003978 -0.0047065 NA 0.0001593 0 NA NA 10 2596.258 -5172.427 45.55036 1.282580e-10
25 -2.836876 -6.9e-06 0.0005439 0 -0.0003425 0.0004039 -0.0013918 2.6e-06 NA NA 0.0004525 -0.0003056 11 2593.413 -5164.719 53.25808 2.718774e-12
19 -2.835390 -6.9e-06 0.0004878 0 -0.0003117 0.0003459 -0.0028671 2.2e-06 0.0000611 NA NA -0.0002915 11 2593.127 -5164.146 53.83075 2.041833e-12
11 -2.836632 -7.1e-06 0.0005408 0 -0.0003235 0.0004522 -0.0042024 1.1e-06 0.0001223 NA 0.0004536 NA 11 2593.075 -5164.042 53.93520 1.937935e-12
5 -2.835167 -7.0e-06 0.0004849 0 -0.0003063 0.0004032 -0.0006555 1.5e-06 NA 0 NA NA 10 2591.819 -5163.548 54.42891 1.514015e-12
28 -2.837380 -7.4e-06 0.0005686 0 -0.0003688 0.0004337 -0.0004705 NA NA 0 0.0004502 -0.0002091 11 2590.242 -5158.376 59.60107 1.140276e-13
14 -2.836906 -7.3e-06 0.0005533 0 -0.0003357 0.0004558 -0.0045421 NA 0.0001503 0 0.0004513 NA 11 2589.975 -5157.843 60.13398 8.735578e-14
22 -2.835720 -7.3e-06 0.0005048 0 -0.0003240 0.0003677 -0.0041769 NA 0.0001335 0 NA -0.0001950 11 2589.900 -5157.692 60.28495 8.100421e-14
27 -2.836861 -7.0e-06 0.0005441 0 -0.0003386 0.0004046 -0.0026825 2.3e-06 0.0000498 NA 0.0004516 -0.0002839 12 2586.845 -5149.562 68.41474 1.390423e-15
21 -2.835393 -6.9e-06 0.0004868 0 -0.0003158 0.0003443 -0.0013506 2.5e-06 NA 0 NA -0.0003208 11 2585.575 -5149.042 68.93502 1.071937e-15
13 -2.836661 -7.1e-06 0.0005419 0 -0.0003330 0.0004604 -0.0007883 1.6e-06 NA 0 0.0004564 NA 11 2585.547 -5148.986 68.99156 1.042059e-15
7 -2.835165 -7.1e-06 0.0004852 0 -0.0002970 0.0003947 -0.0043083 1.0e-06 0.0001339 0 NA NA 11 2585.227 -5148.347 69.62984 7.573414e-16
30 -2.837203 -7.4e-06 0.0005618 0 -0.0003516 0.0004275 -0.0040513 NA 0.0001264 0 0.0004485 -0.0001811 12 2583.611 -5143.095 74.88169 5.481093e-17
29 -2.836865 -6.9e-06 0.0005432 0 -0.0003419 0.0004034 -0.0014547 2.6e-06 NA 0 0.0004525 -0.0003080 12 2579.294 -5134.462 83.51548 7.312432e-19
23 -2.835374 -6.9e-06 0.0004867 0 -0.0003108 0.0003452 -0.0029967 2.2e-06 0.0000624 0 NA -0.0002946 12 2579.012 -5133.896 84.08130 5.510563e-19
15 -2.836649 -7.1e-06 0.0005418 0 -0.0003245 0.0004524 -0.0040732 1.1e-06 0.0001205 0 0.0004536 NA 12 2578.949 -5133.771 84.20563 5.178431e-19
31 -2.836846 -7.0e-06 0.0005430 0 -0.0003377 0.0004039 -0.0027978 2.3e-06 0.0000509 0 0.0004516 -0.0002867 13 2572.729 -5119.310 98.66709 3.749143e-22
kable(Subset.dredged.extinct.model, "html", caption = "BAMM Extinction Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.835512 -7.3e-06 0.0004994 0 -0.0003236 0.0004112 -0.0003583 NA NA NA NA NA 8 2617.018 -5217.977 0.00000 9.982022e-01
8 -2.837029 -7.4e-06 0.0005575 0 -0.0003517 0.0004687 -0.0004496 NA NA NA 0.0004539 NA 9 2610.740 -5203.407 14.56967 6.846268e-04
16 -2.835892 -7.3e-06 0.0005110 0 -0.0003416 0.0003737 -0.0004570 NA NA NA NA -0.0002267 9 2610.652 -5203.231 14.74603 6.268432e-04
2 -2.835374 -7.2e-06 0.0004942 0 -0.0003057 0.0003976 -0.0048344 NA 0.0001612 NA NA NA 9 2610.384 -5202.694 15.28266 4.793257e-04
1 -2.835131 -7.0e-06 0.0004827 0 -0.0003044 0.0004030 -0.0008253 1.5e-06 NA NA NA NA 9 2605.952 -5193.830 24.14720 5.697974e-06
24 -2.837371 -7.4e-06 0.0005678 0 -0.0003683 0.0004333 -0.0005410 NA NA NA 0.0004501 -0.0002114 10 2604.367 -5188.644 29.33283 4.262599e-07
10 -2.836889 -7.3e-06 0.0005522 0 -0.0003347 0.0004556 -0.0046804 NA 0.0001523 NA 0.0004513 NA 10 2604.101 -5188.112 29.86480 3.267087e-07
18 -2.835724 -7.3e-06 0.0005051 0 -0.0003243 0.0003679 -0.0041362 NA 0.0001330 NA NA -0.0001943 10 2604.018 -5187.947 30.03002 3.008033e-07
4 -2.835551 -7.3e-06 0.0005019 0 -0.0003257 0.0004112 -0.0001576 NA NA 0 NA NA 9 2602.884 -5187.695 30.28165 2.652418e-07
17 -2.835405 -6.9e-06 0.0004876 0 -0.0003165 0.0003449 -0.0012820 2.5e-06 NA NA NA -0.0003181 10 2599.694 -5179.298 38.67931 3.982116e-09
9 -2.836626 -7.1e-06 0.0005397 0 -0.0003311 0.0004602 -0.0009541 1.6e-06 NA NA 0.0004565 NA 10 2599.679 -5179.269 38.70807 3.925262e-09
3 -2.835149 -7.1e-06 0.0004842 0 -0.0002961 0.0003945 -0.0044284 1.0e-06 0.0001357 NA NA NA 10 2599.353 -5178.616 39.36092 2.832071e-09
26 -2.837205 -7.4e-06 0.0005620 0 -0.0003517 0.0004276 -0.0040318 NA 0.0001262 NA 0.0004485 -0.0001808 11 2597.730 -5173.352 44.62511 2.037042e-10
12 -2.837067 -7.4e-06 0.0005600 0 -0.0003538 0.0004688 -0.0002498 NA NA 0 0.0004539 NA 10 2596.607 -5173.124 44.85298 1.817677e-10
20 -2.835901 -7.3e-06 0.0005117 0 -0.0003421 0.0003741 -0.0003955 NA NA 0 NA -0.0002247 10 2596.527 -5172.964 45.01281 1.678074e-10
6 -2.835390 -7.3e-06 0.0004952 0 -0.0003067 0.0003978 -0.0047065 NA 0.0001593 0 NA NA 10 2596.258 -5172.427 45.55036 1.282580e-10
25 -2.836876 -6.9e-06 0.0005439 0 -0.0003425 0.0004039 -0.0013918 2.6e-06 NA NA 0.0004525 -0.0003056 11 2593.413 -5164.719 53.25808 2.718774e-12
19 -2.835390 -6.9e-06 0.0004878 0 -0.0003117 0.0003459 -0.0028671 2.2e-06 0.0000611 NA NA -0.0002915 11 2593.127 -5164.146 53.83075 2.041833e-12
11 -2.836632 -7.1e-06 0.0005408 0 -0.0003235 0.0004522 -0.0042024 1.1e-06 0.0001223 NA 0.0004536 NA 11 2593.075 -5164.042 53.93520 1.937935e-12
5 -2.835167 -7.0e-06 0.0004849 0 -0.0003063 0.0004032 -0.0006555 1.5e-06 NA 0 NA NA 10 2591.819 -5163.548 54.42891 1.514015e-12
28 -2.837380 -7.4e-06 0.0005686 0 -0.0003688 0.0004337 -0.0004705 NA NA 0 0.0004502 -0.0002091 11 2590.242 -5158.376 59.60107 1.140276e-13
14 -2.836906 -7.3e-06 0.0005533 0 -0.0003357 0.0004558 -0.0045421 NA 0.0001503 0 0.0004513 NA 11 2589.975 -5157.843 60.13398 8.735578e-14
22 -2.835720 -7.3e-06 0.0005048 0 -0.0003240 0.0003677 -0.0041769 NA 0.0001335 0 NA -0.0001950 11 2589.900 -5157.692 60.28495 8.100421e-14
27 -2.836861 -7.0e-06 0.0005441 0 -0.0003386 0.0004046 -0.0026825 2.3e-06 0.0000498 NA 0.0004516 -0.0002839 12 2586.845 -5149.562 68.41474 1.390423e-15
21 -2.835393 -6.9e-06 0.0004868 0 -0.0003158 0.0003443 -0.0013506 2.5e-06 NA 0 NA -0.0003208 11 2585.575 -5149.042 68.93502 1.071937e-15
13 -2.836661 -7.1e-06 0.0005419 0 -0.0003330 0.0004604 -0.0007883 1.6e-06 NA 0 0.0004564 NA 11 2585.547 -5148.986 68.99156 1.042059e-15
7 -2.835165 -7.1e-06 0.0004852 0 -0.0002970 0.0003947 -0.0043083 1.0e-06 0.0001339 0 NA NA 11 2585.227 -5148.347 69.62984 7.573414e-16
30 -2.837203 -7.4e-06 0.0005618 0 -0.0003516 0.0004275 -0.0040513 NA 0.0001264 0 0.0004485 -0.0001811 12 2583.611 -5143.095 74.88169 5.481093e-17
29 -2.836865 -6.9e-06 0.0005432 0 -0.0003419 0.0004034 -0.0014547 2.6e-06 NA 0 0.0004525 -0.0003080 12 2579.294 -5134.462 83.51548 7.312432e-19
23 -2.835374 -6.9e-06 0.0004867 0 -0.0003108 0.0003452 -0.0029967 2.2e-06 0.0000624 0 NA -0.0002946 12 2579.012 -5133.896 84.08130 5.510563e-19
15 -2.836649 -7.1e-06 0.0005418 0 -0.0003245 0.0004524 -0.0040732 1.1e-06 0.0001205 0 0.0004536 NA 12 2578.949 -5133.771 84.20563 5.178431e-19
31 -2.836846 -7.0e-06 0.0005430 0 -0.0003377 0.0004039 -0.0027978 2.3e-06 0.0000509 0 0.0004516 -0.0002867 13 2572.729 -5119.310 98.66709 3.749143e-22
#Run model for DR
Subset.MCC.DR.top <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.DR.top, 'data/Subset.MCC.DR.top.rds')

#Run model for ND
Subset.MCC.ND.top <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.ND.top, 'data/Subset.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Subset.data.noMCC <- SS.subset %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP)

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Subset.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.DR.model.data <- lapply(Subset.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Subset.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(SS.subset$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Subset.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(0.9711, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.DR.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.DR.pgls.models, "data/Subset.DR.pgls.models.rds")

#Now for Node Density:
Subset.ND.model.data <- lapply(nd.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.ND.model.data <- lapply(Subset.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.ND.model.data, y = Subset.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.ND.pgls.models, "data/Subset.ND.pgls.models.rds")

#BAMM Top Models 

Subset.MCC.Lambda.top <- gls(log(mean.lambda) ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Lambda.top, 'data/Subset.MCC.Lambda.top.rds')

Subset.MCC.Mu.top <- gls(log(mean.mu) ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Mu.top, 'data/Subset.MCC.Mu.top.rds')

#Now for BAMM 100 models

Subset.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.BAMM.model.data <- lapply(Subset.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.lambda.pgls.models, "data/Subset.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Subset.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ Sexual_selection_ppca
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.mu.pgls.models, "data/Subset.BAMM.mu.pgls.models.rds")

The phylogenetic signal from the DR model is ~0.97. The other models had phylogenetic signals ~ 1, as such they were fixed at 1 to avoid problems of convergence.

Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')
Subset.MCC.DR.top[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.9711
Subset.DR.pgls.models <- readRDS("data/Subset.DR.pgls.models.rds")
Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')

Subset.cols <- brewer.pal(n = 7, name = "Dark2")[-6]

Subset.DR.pgls.summary <- lapply(Subset.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.DR.summary <- data.frame(Subset.MCC.DR.top$coefficients, confint(Subset.MCC.DR.top)) %>% tibble::rownames_to_column()

Subset.DR.pgls.summary <- bind_rows(Subset.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Sexual_selection_ppca` = "Sexual Selection"
                    )

Subset.DR.pgls.summary$Parameter = factor(Subset.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.DR.summary$Parameter = factor(Subset.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.DR.plot <-Subset.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.DR.plot <- Subset.DR.plot + geom_errorbarh(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS DR")

####################################
#For ND

Subset.ND.pgls.models <- readRDS("data/Subset.ND.pgls.models.rds")
Subset.MCC.ND.top <- readRDS('data/Subset.MCC.ND.top.rds')

Subset.ND.pgls.summary <- lapply(Subset.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.ND.summary <- data.frame(Subset.MCC.ND.top$coefficients, confint(Subset.MCC.ND.top)) %>% tibble::rownames_to_column()

Subset.ND.pgls.summary <- bind_rows(Subset.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.ND.pgls.summary$Parameter = factor(Subset.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.ND.summary$Parameter = factor(Subset.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.ND.plot <-Subset.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.ND.plot <- Subset.ND.plot + geom_errorbarh(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS ND")

######################

#For Lambda
Subset.BAMM.lambda.pgls.models <- readRDS("data/Subset.BAMM.lambda.pgls.models.rds")
Subset.MCC.Lambda.top <- readRDS('data/Subset.MCC.Lambda.top.rds')

Subset.Lambda.pgls.summary <- lapply(Subset.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Lambda.summary <- data.frame(Subset.MCC.Lambda.top$coefficients, confint(Subset.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Subset.Lambda.pgls.summary <- bind_rows(Subset.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Lambda.pgls.summary$Parameter = factor(Subset.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Lambda.summary$Parameter = factor(Subset.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Lambda.plot <-Subset.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Lambda.plot <- Subset.Lambda.plot + geom_errorbarh(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Speciaton")

######################

#For Mu

Subset.BAMM.Mu.pgls.models <- readRDS("data/Subset.BAMM.mu.pgls.models.rds")
Subset.MCC.Mu.top <- readRDS('data/Subset.MCC.Mu.top.rds')

Subset.Mu.pgls.summary <- lapply(Subset.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Mu.summary <- data.frame(Subset.MCC.Mu.top$coefficients, confint(Subset.MCC.Mu.top)) %>% tibble::rownames_to_column()

Subset.Mu.pgls.summary <- bind_rows(Subset.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Mu.pgls.summary$Parameter = factor(Subset.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Mu.summary$Parameter = factor(Subset.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Mu.plot <-Subset.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Mu.plot <- Subset.Mu.plot + geom_errorbarh(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Extinction")

#########################################
grid.arrange(symmetrise_scale(Subset.DR.plot, "x"),
             symmetrise_scale(Subset.ND.plot, "x"),
             symmetrise_scale(Subset.Lambda.plot, "x"), 
             symmetrise_scale(Subset.Mu.plot, "x"), 
             nrow = 1)

Figure S11: Measures of sexual selection a phylogenetic pca of sexual dimorphism, social polygyny and paternal care provide some evidence to suggest sexual selection predicts extinction rate (DR and ND models). This figure is seen in Figure 1 within the manuscript

Table S11: Highest posterior density intervals are calculated for the above figure (ie; models using Male-bias sexual selection measures). using this measure the models using \(\lambda_{DR}\) have model estimates that fall on the positive side 95 % of the time in 100 trees. For \(\lambda_{ND}\), there is a positive skew, however the 95 % HPD interval overlaps zero. \(\lambda_{BAMM}\) also shows a lesser positive skew. The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7 and Table S9.

Subset.hpd.DR.top <- list()
Subset.DR.pgls.summary <- na.omit(Subset.DR.pgls.summary)
for(x in unique(Subset.DR.pgls.summary$Parameter)) {
Subset.hpd.DR.top[[x]] = hdi(Subset.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.DR.top <- bind_rows(Subset.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Subset.hpd.DR.top, 'data/Subset.hpd.DR.top.rds')

#For ND
Subset.hpd.ND.top <- list()
Subset.ND.pgls.summary <- na.omit(Subset.ND.pgls.summary)
for(x in unique(Subset.ND.pgls.summary$Parameter)) {
Subset.hpd.ND.top[[x]] = hdi(Subset.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.ND.top <- bind_rows(Subset.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.ND.top, 'data/Subset.hpd.ND.top.rds')

###############
Subset.hpd.Lambda.top <- list()
Subset.Lambda.pgls.summary <- na.omit(Subset.Lambda.pgls.summary)
for(x in unique(Subset.Lambda.pgls.summary$Parameter)) {
Subset.hpd.Lambda.top[[x]] = hdi(Subset.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Lambda.top <- bind_rows(Subset.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Subset.hpd.Lambda.top, 'data/Subset.hpd.Lambda.top.rds')

################
Subset.hpd.Mu.top <- list()
Subset.Mu.pgls.summary <- na.omit(Subset.Mu.pgls.summary)
for(x in unique(Subset.Mu.pgls.summary$Parameter)) {
Subset.hpd.Mu.top[[x]] = hdi(Subset.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Mu.top <- bind_rows(Subset.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.Mu.top, 'data/Subset.hpd.Mu.top.rds')

Subset.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Subset DR HPD Interval")
Subset DR HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower 0.00911 -0.0158 -9.31e-05 0.000313 -0.0036 -2.96e-06
Upper 0.0608 -0.000501 2.76e-05 0.0128 0.00894 1.19e-06
Subset.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Subset ND HPD Interval")
Subset ND HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00031 -0.000426 -3.63e-06 -3.05e-05 -0.000249 -9.23e-08
Upper 0.00154 0.000166 1.64e-06 0.000487 0.000249 6.61e-08
Subset.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Speciation HPD Interval")
Subset BAMM Speciation HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0198 -0.00853 -7.27e-05 -0.0114 -0.0106 -1.23e-06
Upper 0.0313 0.00801 5.36e-05 0.00827 0.00395 1.91e-06
Subset.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Extinction HPD Interval")
Subset BAMM Extinction HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0441 -0.0134 -0.000145 -0.0125 -0.0156 -1.97e-06
Upper 0.0431 0.013 7.66e-05 0.0104 0.00951 3.36e-06

Additional Figures

MCC.BAMM.model.data <- readRDS('data/MCC.BAMM.model.data.rds')
restricted.data %>% ggplot(aes(x = MCC.DR, y = MCC.BAMM.model.data$mean.lambda))+
  geom_point()+
  scale_y_continuous(trans = "log")+
  theme_minimal()

R Session information

sessionInfo() %>% pander

R version 3.3.1 (2016-06-21)

**Platform:** x86_64-apple-darwin13.4.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), ggtree(v.1.6.11), EBImage(v.4.16.0), RColorBrewer(v.1.1-2), kableExtra(v.0.7.0), phangorn(v.2.4.0), HDInterval(v.0.2.0), MuMIn(v.1.42.1), coda(v.0.19-1), car(v.3.0-2), carData(v.3.0-1), BAMMtools(v.2.1.6), dplyr(v.0.7.6), ggExtra(v.0.8), reshape2(v.1.4.3), purrr(v.0.2.4), caper(v.1.0.1), mvtnorm(v.1.0-6), MASS(v.7.3-50), ggridges(v.0.5.0), brms(v.2.4.0), Rcpp(v.0.12.18), phytools(v.0.6-44), maps(v.3.2.0), gridExtra(v.2.3), geiger(v.2.0.6), lme4(v.1.1-18-1), Matrix(v.1.2-14), repmis(v.0.5), diversitree(v.0.9-10), ape(v.5.1), mgcv(v.1.8-24), nlme(v.3.1-128), tidyr(v.0.7.2), knitr(v.1.20), pander(v.0.6.2) and ggplot2(v.3.0.0)

loaded via a namespace (and not attached): readxl(v.1.0.0), backports(v.1.1.2), fastmatch(v.1.1-0), plyr(v.1.8.4), igraph(v.1.1.2), lazyeval(v.0.2.1), splines(v.3.3.1), crosstalk(v.1.0.0), rstantools(v.1.5.1), inline(v.0.3.15), digest(v.0.6.16), htmltools(v.0.3.6), tiff(v.0.1-5), rsconnect(v.0.8.8), gdata(v.2.18.0), magrittr(v.1.5), openxlsx(v.4.1.0), readr(v.1.1.1), matrixStats(v.0.54.0), R.utils(v.2.7.0), xts(v.0.11-0), jpeg(v.0.1-8), colorspace(v.1.3-2), rvest(v.0.3.2), haven(v.1.1.2), jsonlite(v.1.5), bindr(v.0.1.1), survival(v.2.39-4), zoo(v.1.8-3), glue(v.1.3.0), gtable(v.0.2.0), R.cache(v.0.13.0), rstan(v.2.17.3), BiocGenerics(v.0.20.0), abind(v.1.4-5), scales(v.1.0.0), miniUI(v.0.1.1.1), plotrix(v.3.7-3), viridisLite(v.0.3.0), xtable(v.1.8-2), foreign(v.0.8-66), subplex(v.1.4-1), deSolve(v.1.20), stats4(v.3.3.1), StanHeaders(v.2.17.2), DT(v.0.4), animation(v.2.5), htmlwidgets(v.1.2), httr(v.1.3.1), threejs(v.0.3.1), gplots(v.3.0.1), pkgconfig(v.2.0.2), loo(v.2.0.0), R.methodsS3(v.1.7.1), locfit(v.1.5-9.1), labeling(v.0.3), tidyselect(v.0.2.4), rlang(v.0.2.2), later(v.0.7.3), munsell(v.0.5.0), cellranger(v.1.1.0), tools(v.3.3.1), fftwtools(v.0.9-8), evaluate(v.0.10.1), stringr(v.1.2.0), yaml(v.2.2.0), zip(v.1.0.0), caTools(v.1.17.1), mime(v.0.5), R.oo(v.1.22.0), xml2(v.1.2.0), bayesplot(v.1.6.0), shinythemes(v.1.1.1), png(v.0.1-7), curl(v.3.2), clusterGeneration(v.1.3.4), tibble(v.1.3.4), stringi(v.1.1.6), Brobdingnag(v.1.2-6), forcats(v.0.3.0), lattice(v.0.20-33), nloptr(v.1.0.4), markdown(v.0.8), shinyjs(v.1.0), msm(v.1.6.6), combinat(v.0.0-8), bridgesampling(v.0.5-2), data.table(v.1.11.4), bitops(v.1.0-6), httpuv(v.1.4.5), R6(v.2.2.2), promises(v.1.0.1), KernSmooth(v.2.23-15), rio(v.0.5.10), codetools(v.0.2-14), colourpicker(v.1.0), gtools(v.3.8.1), assertthat(v.0.2.0), rprojroot(v.1.3-2), withr(v.2.1.2), shinystan(v.2.4.0), mnormt(v.1.5-5), expm(v.0.999-2), hms(v.0.4.2), quadprog(v.1.5-5), minqa(v.1.2.4), rmarkdown(v.1.10), numDeriv(v.2016.8-1), scatterplot3d(v.0.3-41), shiny(v.1.1.0), base64enc(v.0.1-3) and dygraphs(v.1.1.1.6)

References

Armenta, J. K., P. O. Dunn, and L. A. Whittingham. 2008. Quantifying avian sexual dichromatism: A comparison of methods. Journal of Experimental Biology 211:2423. Journal Article.

BirdLife International and Handbook of the Birds of the World. 2017. Bird species distribution maps of the world. http://datazone.birdlife.org/species/requestdis.

Dale, J., C. J. Dey, K. Delhey, B. Kempenaers, and M. Valcu. 2015. The effects of life history and sexual selection on male and female plumage colouration. Nature 527:367–370. Journal Article.

Del Hoyo, J., A. Elliott, and D. Christie. n.d. Handbook of the birds of the world (Vols. 8-16). Book, Lynx Edicions 2003-2011.

Harvey, M. G., G. F. Seeholzer, B. T. Smith, D. L. Rabosky, A. M. Cuervo, and R. T. Brumfield. 2017. Positive association between population genetic differentiation and speciation rates in new world birds. Proceedings of the National Academy of Sciences 114:6328–6333.

Harvey Michael, G., L. Rabosky Daniel, and N. Cooper. 2017. Continuous traits and speciation rates: Alternatives to state‐dependent diversification models. Methods in Ecology and Evolution 9:984–993. Journal Article.

Jetz, W., G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. 2012. The global diversity of birds in space and time. Nature 491:444–448. Journal Article.

Quintero, I., and W. Jetz. 2018. Global elevational diversity and diversification of birds. Nature 555:246.

Schliep, K. P. 2010. Phangorn: Phylogenetic analysis in r. Bioinformatics 27:592–593.